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EXACT  PERFORMANCE  OF  FILTERED  AND  WEIGHTED  ENERGY  DETECTOR  WITH 
MISMATCHED  FREQUENCY  AND  TIME  LOCATIONS  AND  CHARACTERISTICS 

INTRODUCTION 

The  detection  of  weak  random  signals  of  unknown  duration, 
bandwidth,  time  location,  and  frequency  shift  in  the  presence  of 
colored  noise  is  often  accomplished  in  practice  by  filtering  the 
received  waveform  to  the  frequency  band  of  interest,  squaring  the 
filter  output,  and  time-weighting  this  quantity  prior  to 
accumulation  and  threshold  comparison.  The  performance  of  this 
processor  obviously  depends  on  the  spectra  of  the  received  signal 
and  noise  processes  as  well  as  on  the  transfer  characteristics  of 
the  filter  employed.  Lack  of  detailed  knowledge  of  the  signal 
spectral  characteristics  or  frequency  shift  will  often  dictate 
that  a  fairly  broad  (mismatched)  filter  passband  be  utilized  in 
order  that  the  signal  be  passed  when  present.  Similarly, 
uncertainty  about  the  signal  location  or  duration  will  mandate 
that  the  receiver  observation  time  be  lengthened  in  order  to 
ensure  capture  of  any  impingent  signal  energy. 

Additionally,  the  received  signal  process  is  often 
characterized  as  a  finite-duration  burst  of  a  stationary  process. 
Since  the  signal  time  location  and  duration  are  often  unknown, 
the  receiving  filter  must  be  kept  open  at  all  times,  thereby 
allowing  the  filter  output  noise  to  reach  its  full-strength 
steady  state  value.  By  contrast,  a  gated  input  signal  leads  to  a 
filter  output  component  which  gradually  builds  up  in  time  and 
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then  decays  slowly  after  the  signal  is  turned  off.  These 
nonstationary  effects,  coupled  with  the  particular  time  weighting 
employed  during  the  observation  interval,  influence  the 
performance  of  the  detector  in  a  complicated  nonobvious  fashion. 

There  is  a  need  to  be  able  to  calculate  exactly  the 
performance  of  this  type  of  signal  processor  in  the  presence  of 
such  deleterious  factors,  so  that  the  degradations  associated 
with  lack  of  knowledge  can  be  ascertained  quantitatively. 
Furthermore,  it  will  not  suffice  to  resort  to  Gaussian 
approximations  for  the  processor  output,  because  the  number  of 
samples  employed  are  not  large  enough  to  utilize  the  central 
limit  theorem,  especially  on  the  tails  of  the  distribution  (small 
false  alarm  probabilities)  where  we  are  interested.  This  need  is 
addressed  in  this  report  for  the  case  of  Gaussian  input  signals 
in  Gaussian  noise,  with  the  result  that  programs  are  furnished 
for  exactly  assessing  the  performance  of  the  mismatched  energy 
detector  for  a  very  wide  variety  of  characteristics  and 
selections  of  parameters.  In  particular,  no  presumptions  are 
made  about  the  input  signal  time-bandwidth  product  or  about  the 
size  of  the  product  of  observation  time  and  filter  bandwidth. 
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SYSTEM  DESCRIPTION  AND  ASSUMPTIONS 


Figure  1.  Block  Diagram  of  Detector 


The  detector  of  interest  is  the  digital  processor  depicted  in 
figure  1;  time  sampling  increment  A  is  arbitrary  but  should  be 
approximately  matched  to  the  coherence  times  of  the  input  signal 
and  filter.  The  digital  input  sequence  x(kA),  filter  impulse 
response  A  h(mA),  and  summer  weights  w{kA)  are  all  real.  Input 
x(kA)  consists  of  either  noise-alone  or  gated  signal  plus  noise, 
according  to  model 


X  { kA ) 


n(  kA ) 
or 

is(kA)  g(kA)  +  n(kA) 


>  for  all  k 


(1) 


Input  noise  n  kA )  is  present  for  all  time  and  is  a  stationary 
zero-mean  discrete  Gaussian  sequence  with  covariance 

R^(mA)  =  n(kA)  n(kA-mA)  for  all  m,k  ,  (2) 

where  an  overbar  denotes  an  ensemble  average.  The  numerical 
evaluation  of  the  values  of  the  noise  covariance,  directly  from 
a  specified  noise  spectrum  Gj^(f),  is  considered  in  appendix  A. 
Underlying  input  signal  s{kA)  in  (1)  (if  present)  is  also  a 
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stationary  zero-mean  discrete  v^aussian  sequence  with  covariance 

Rg{mA)  =  s(kA)  s(kA-raA)  for  all  m,k  .  (3) 

However,  input  signal  samples  s(kA)  are  gated  by  function  g{kA) 
which  is  nonzero  only  in  a  limited  time  region: 

g(kA)  =  1  for  k^  i  k  i  kj^  ,  zero  otherwise  .  (4) 

This  results  in  a  gated  burst  of  stationary  signal  sequence  s(kA) 
being  inputted  to  digital  filter  h{mA);  the  input  signal 
starting  and  ending  times  k  A  and  k.  A,  respectively,  are 

3.  O 

generally  unknown,  except  in  an  approximate  way.  This  generality 
allows  for  consideration  of  input  signals  of  unknown  arrival  time 
and  duration  at  the  detector  input. 

The  filter  output  y(kA)  can  be  conveniently  broken  into 
signal  and  noise  components,  in  accordance  with  (1),  and  is 
available  by  means  of  discrete  convolution 

y(kA)  =  X]  A  h(mA)  x(kA  -  mA)  =  yg(kA)  +  y^^CkA)  ,  (5) 

m 

when  an  input  signal  is  present.  (Summations  without  limits  are 
over  (-“>,+“>).)  Howeve’',  due  to  the  gating  in  (1),  filter  output 
signal  yg(kA)  will  have  transient  phases,  including  a  buildup 
just  after  time  k^A  and  a  decay  just  after  time  k^^A.  This 
nonstationary  behavior  is  included  and  exactly  accounted  for  in 
the  following  analysis. 

The  filter  output  (5)  is  then  squared,  scaled  by  weights 
w(kA).  and  summed  to  give  system  output 
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z  =  H  w(kA)  y^(kA)  ,  (6) 

k 

where  the  weights  are  nonzero  only  in  a  limited  observation  time 
specified  by 

w(kA)  ^  0  for  k^  i  k  <:  k^  .  (7) 

These  weights  need  not  be  uniform;  for  example,  they  could  be 

exponential.  Output  z  is  compared  with  a  threshold  for  a 

declaration  of  signal  absent  versus  signal  present. 

Observation  time  limits  k  A  and  k,A  are  under  the  control  of 

c  d 

the  receiver  processor,  and  would  hopefully  match  fairly  well  the 
time  interval  over  which  the  gated  input  signal  is  received;  see 
(4).  But,  in  any  event,  the  degree  of  generality  in  (7)  allows 
for  investigation  of  observation  times  that  might  be  too  short, 
thereby  losing  some  signal  energy,  or  too  long,  thereby  picking 
up  additional  unwanted  noise  components;  both  situations  lead  to 
deterioration  of  the  detectability  of  weak  signals  and  should  be 
avoided  if  possible.  An  illustration  of  the  parameters  is  given 
in  figure  2,  for  signal-only  present. 


input  signal  samples 
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Figure  2.  Gated  and  Observation  Intervals;  No  Noise 
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FILTER  OUTPUT  NOISE  CHARACTERIZATION 

The  filter  output  noise  sequence  is  obtained  from  discrete 
convolution  (5)  and  (1)  as 

y  (kA)  =  A  h(mA)  n(kA  -  mA)  for  all  k  ,  (8) 

m 

where  it  is  presumed  that  the  filter  h(mA)  has  been  open  to  noise 
input  n(kA)  for  all  time;  this  results  in  noise  output  yj^(kA) 
being  a  zero-mean  stationary  Gaussian  sequence.  The  filter 
output  noise  covariance  at  general  times  kj^A  and  k2A  is  given  by 

Cn(kj/S,kjiS)  -  y„(k,A)  y„(k2fi)  = 

=  )  h(mfi)  h(pdj  n(kjd  -  md)  n(k2fl  -  pfl)  = 

mp 

=  YZ.  h(mA)  h(pA)  Rj^(kjA  -  k2A  +  pA  -  mA )  = 
mp 


=  II  ^n^’^l^  -  ^2^  ■ 

where  the  filter  correlation  function  is  defined  as 

(^j^(jA)  =  tZ  I  h(mA)  h{mA  -  jA)  for  all  j  .  (10) 

m 

The  right-hand  side  of  (9)  is  a  function  only  of  the  time 
difference  k^^A  -  k2A,  in  keeping  with  the  stationarity  of  filter 
output  noise  y^(kA). 
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If  the  input  noise  in  figure  1  is  white,  its  covariance  in 
( 2 )  becomes 


'mA ) 


(11) 


2 

where  is  the  input  noise  power.  In  this  case,  the  filter 
output  noise  covariance  in  (9)  simplifies  to 


Cj^(kj^A,k2A)  =  +j^(kj^A  -  k2A)  for  white  input  noise 


(12) 


The  numerical  evaluation  of  filter  correlation  function  ♦j^(jA)  in 
(10),  directly  from  a  specified  transfer  function  H(f),  is 
considered  in  appendix  B.  A  possible  problem  with  discontinuous 
functions  is  treated  in  appendix  C. 
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FILTER  OUTPUT  SIGNAL  CHARACTERIZATION 

The  filter  output  signal  sequence  (when  present)  follows  from 
( 5 )  and  ( 1 )  as 

yg(kA)  =  A  h(kA  -  mA)  s(mA)  g(mA)  for  all  k  .  (13) 

m 

The  finite  duration  of  gating  sequence  g(kA)  in  (4),  as  well  as 
the  realizability  of  impulse  response  A  h{mA),  will  serve  to 
terminate  the  summation  in  (13)  at  finite  limits. 

We  presume  that  filter  h(mA)  in  figure  1  is  realizable;  that 
is, 

h(mA)  =  0  for  m  <  0  .  (14) 

This  makes  filter  output  signal  y_(kA)  in  (13)  equal  to  zero  for 

k  <  k  .  On  the  other  hand,  if  the  filter  has  an  infinite- 
a 

duration  impulse  response  A  h(mA),  then  y  (kA)  is  nonzero  for 
k  i  k^;  see  figure  3.  However,  the  filter  output  noise  yj^(kA) 
will  dominate  the  signal  output  for  k  >>  k^^,  meaning  that  per¬ 
formance  of  the  detector  will  be  poor  in  the  case  when  k^  >>  k^^. 
That  is,  too  long  an  observation  time  in  (7),  relative  to  the 
signal  duration,  is  detrimental  to  signal  detectability. 

steady  state 


buildup  91  T  decay 
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The  signal  output,  yg(kA)  in  (13),  from  filter  h(mA)  is  not 
stationary.  The  filter  output  signal  covariance  is  defined  as 


Cg(k^A,k2A)  =  yg(k^A)  yg(k2A)  =  Cg{k2A,k^A)  for  all  k^,k2.  (15) 

This  function  is  zero  for  kj^  <  k^  or  k2  ^  therefore,  in  the 

following,  we  can  confine  calculation  of  (15)  to  kj^  1  k^  and 
k„  1  k^.  The  particular  filter  output  signal  sample  y^(k^A)  is 

Z  CL  S  ol 

nonzero  only  if  h(0)  ^  0. 

Substitution  of  (13)  in  (15)  and  the  use  of  (3)  yields 


C  (kj^A,k2A)  =  A"^  )  [  h(kjA-mA)  h(k2A-pA)  s(mA)  s(pA)  g(mA)  g(pA)  = 

mp 


=  A  )  h(k,A-mA)  h(k-A-pA)  R  (mA-pA)  g(mA)  g(pA)  .  (16) 

mp 


When  we  take  explicit  account  of  realizability  condition  (14)  and 
the  finite  duration  of  unity  gating  function  g(kA)  in  (4),  the 
filter  output  signal  covariance  in  (16)  can  be  efficiently 
computed  according  to 


K. 


K. 


Cg(kjA,k2A)  =  A^  ^  h(kjA-mA)  ^  h(k2A-pA)  Rg(mA-pA)  ,  (17) 


m=k 


p=k 


where 


=  min(kj,kjj)  ,  K2  =  min(k2,kj^) 


(18) 


Since  k.  i  k^,  k,  i  k,,  k-  i  k^,  it  follows  that  K,  i.  k,  and 

OoiaZa  la 

K-  i.  k^,  thereby  making  the  summations  in  (17)  treat  only 

Z  a 

nonzero  entries. 
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SUMMER  CONSIDERATIONS 

The  weighted  sununer  in  figure  1  and  (6)  -  (7)  takes  its 

samples  at  times  kA,  where  k^  i  k  i  k^.  If  k^  <  k^,  then 

signal  gets  into  the  summer ,  and  the  signal  at  the  input  cannot 

be  detected  at  all,  whether  present  or  not.  Since  this  is  not  a 

useful  application  of  the  detector  in  figure  1,  the  only 

situation  we  will  address  is  that  where  kj  1  k_ ,  so  that  some 

a  a  - 

signal  (when  present)  contributes  to  the  summer  output. 

There  are  then  two  possibilities  for  observation  start  time 

k^A,  as  indicated  in  figure  4  below.  If  k^  <  k^  (case  i),  the 

summer  is  taking  in  noise-only  samples  for  k^  1  k  <  k,;  this  will 

degrade  performance  of  the  detector  but  is  sometimes  unavoidable 

when  the  signal  onset  time,  k^^A,  is  unknown.  In  this  case,  we 

only  need  to  compute  signal  output  covariance  Cg(kj^A,k2A)  for 

k  1  k.,kT  1  kj,  since  C„  =  0  for  k,  <  k^  or  k-  <  k,. 
al  2d  s  la  2a 

On  the  other  hand,  if  k  1  k^  (case  2),  signal  is  taken  into 

Q  cl 

the  summer  on  all  samples  available  to  it.  However,  if  k^  is 
considerably  larger  than  k^,  a  significant  portion  of  the  signal 
contribution  can  be  lost;  this  degradation  of  performance  is 


'cl 


—  ygCkA) 


c2 


Figure  4.  Filter  Output  Signal  y  (kA) 


10 


TR  8913 


sometimes  unavoidable  when  k_  is  unknown.  In  this  latter  case, 

a 

we  only  need  to  compute  Cg{kj^A,k2A)  for  k^  i  ^  ^d* 

The  general  rule  is  that  we  need  to  compute,  from  (17),  the 
filter  output  signal  covariance  Cg{k^A,k2A)  for 

s  max(k^,k^)  1  k^,k2  1  k^  .  (19) 

Then,  if  k  <  k^,  use  0  for  C_(k,A,k-5A)  for  those  values  where 

k,  <  k  or  k-,  <  k^.  Of  course,  noise  covariance  C  (k,A,k,A)  in 
X  £1  ^  o  n  X  z 

(9)  or  (12)  must  be  computed  for  k^  1  ^  ^d  cases. 

.‘Jinre,  from  ( 19  )  , 

k.  i.  max(k,,k^)  1  k^,  then  K.  =  min(k.,k.  )  i  k  ;  (20) 

X  oCo  X  XOo 

also,  since 

k2  i  max(k^,k^)  i.  k^,  then  K2  =  min(k2,kjj)  i  k^  .  (21) 

Therefore,  the  sums  in  (17)  always  have  some  entries;  that  is, 
the  upper  limit  is  never  less  than  the  lower  limit.  The  only 
restriction  is 

i  .  (22) 

Of  course,  we  must  always  have  k^  i  k^^  and  k^  i  k^. 

In  summary,  the  filter  output  signal  covariance  follows  from 
(17)  and  (18),  while  the  filter  output  noise  covariance  is 
available  in  (9)  or  (12).  The  region  where  the  filter  output 
signal  covariance  Cg(kj^A,k2A)  calculation  is  nonzero  is 
crosshatched  in  figure  5  for  the  case  where  k^  <  k_.  If 
^  then  C„(k,A,k-,A)  is  nonzero  in  the  larger  square 

C  u  S  X  z 

indicated  in  figure  5. 
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Figure  5.  Calculation  Regie 
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CHARACTERISTIC  FUNCTION  OF  SYSTEM  OUTPUT  z 

The  system  output  z  is  given  in  (6)  as  a  sum  of  weighted  and 
squared  correlated  zero-mean  Gaussian  random  variables  y(kA). 
Also,  the  covariance  of  the  noise  component  yj^(kA)  is  given  by 
(9)  or  (12),  while  the  covariance  of  the  signal  component  yg(kA) 
(if  present)  is  given  by  (17)  -  (18)  and  figure  5.  Therefore, 
the  covariance  of  y(kA)  is  given  by 

“  Cg(kjA,k2A)  +  Cj^ ( k A , k2A )  .  (23) 

The  system  input  signal-to-noise  ratio  is  R  (0)/R  (0)  in  terms  of 

9  11 

input  covariances  (2)  and  (3). 

For  nonnegative  weights  {w(kA)},  we  define  random  variables 

aj^  =  |w ( kA y  y(kA)  for  k^  i  k  1  k^  .  (24) 

Then,  (6)  -  (7)  yields  the  output  in  the  form 

kd 

z  =  YU.  f  (25) 

k=k^ 

where  zero-mean  Gaussian  random  variables  (aj^j  have  covariance 

^a^^l'^2^  =  ®k^^^  ^  ^w(kjA)  w(k2A)  Cy(kj^A,k2A)  .  (26) 

In  order  to  find  the  characteristic  function  of  random 
variable  z  in  (25),  we  consider  the  square  symmetric  covariance 
matrix  C  with  elements  (26)  for  k^  1  ^  ^d*  ^ 
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normalized  modal  matrix  of  the  orthonormal  eigenvectors  of  C  and 
let  A  be  the  diagonal  matrix  of  the  eigenvalues  of  C;  see 

[1;  pages  36  -  39].  Then  we  have 

Q  =  Q  q'^  =  I  /  C  Q  =  A  .  (27) 

Now  let  column  vector  A  be  made  up  of  elements  {aj^j  for 
1  k  i  k^,  as  defined  in  (24),  in  which  case  (25)  and  (26)  can 
be  expressed  as 

z  =  A  ,  C  =  A  a"^.  (28) 

Also,  let  column  vector  B  be  defined  by  linear  transformation 

B  =  A  ;  then  A  =  Q  B  .  (29) 

Substitution  of  (29)  into  (28)  yields 

kd 

z  =  a"^  A  =  B*^  Q  B  =  b"^  B  =  b^  ,  (30) 

k=k  ^ 
c 

where  {bj^j  are  the  elements  of  B.  At  the  same  time,  the 
covariance  matrix  of  vector  B  is 

B  B^  =  q'^  A  a"^  Q  =  C  Q  =  A  ,  (31) 

where  we  used  (29),  (28),  and  (27).  Thus,  the  sum  for  z  in  (30) 
is  composed  of  uncorrelated  (and  therefore  independent)  zero-mean 
Gaussian  random  variables  where  the  variance  of  element  bj^ 

is  eigenvalue  Xj^. 

The  characteristic  function  of  random  variable  z  in  (30)  is. 
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using  the  independence  and  Gaussian  character  of 

=  exp(  i?;zT  =  I"  expjub^  = 

k=k 

c 


-ij 


=  I  du  exp(ii;u‘')  (2nX,^)  ^  exp 


k=k 


-u 


2X, 


Li 

=n( 

k=k 


1  -  i2?;x. 


-h 


(32) 


Although  this  final  result  is  given  as  a  finite  product  of 
k^-k^+1  principal-value  square  roots,  it  can  be  computed  as  a 
single  principal-value  square  root  of  a  finite  product  of  complex 
first-order  polynomials,  provided  that  the  location  of  the 
product  in  the  complex  plane  is  tracked  from  ^  *  0;  for  example, 
see  [2;  pages  B-1  and  B-2]. 
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COVARIANCE  AND  FILTER  EXAMPLES 


INPUT  SIGNAL  CHARACTERIZATION 


The  input  signal  covariance  will  be  taken  to  be  a  sum  of 
damped  exponentials: 


J 

“-i  ^  /  (33) 

s  j=l  :  SD 


where  the  J  exponential  parameters  l^gj)  are  all  distinct.  The 
scalings  {a- j  and  time  constants  (t  • )  can  be  complex,  provided 
that  they  occur  in  complex  conjugate  pairs;  that  is,  input  signal 
covariance  Rg(T)  must  be  real  for  all  t.  The  origin  value. 


J 

R.iO)  =  YU 
s  ^  1 


(34) 


is  the  input  power  of  stationary  signal  sequence  s ( kA )  in  (1)  and 
(3),  prior  to  gating.  The  input  signal  spectrum  corresponding 
to  (33)  is,  with  w  =  2nf, 


Gs(f)  =  I  dT  exp(-i2nfr)  R^{t)  =  ^  ^ ,  (35) 

3=1  1  +  (W  Tgj ) 

which  is  seen  to  contain  2J  distinct  poles  in  the  complex 
f-plane.  ^^(f)  can  contain  up  to  2J-2  zeros. 

For  time  sampling  increment  A,  the  sampled  input  signal 
covariance  is,  from  (33), 
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RgCJtA)  =  YIZ  exp(-a .  |k|  )  ;  a.  =  ^  . 

®  j  =  l  ^  ^  ^  "^sj 


The  J  exponential  parameters  {a^}  are  distinct,  but  they  can  be 
complex;  in  fact,  {Oj}  and  {a^j  can  both  be  complex  provided  that 
they  occur  in  complex  conjugate  pairs  such  that  Rg(kA)  is  real 
for  all  k.  In  fact,  some  of  the  scalings  {Oj)  can  be  negative, 
provided  that  total  spectrum  Gg(f)  in  (35)  is  nonnegative  for  all 
f.  The  exponential  parameters  {a^}  in  (36)  must  also  satisfy 
Re  aj  >0  in  order  that  the  covariance  tend  to  zero  at  k  =  I®. 

Fcr  a  particular  pair  of  complex  conjugate  terms,  say 


a,  =  a  ,  a.,  =  ,  t  ,,  =  t  ,  a,  =  a  ,  a~  =  a  , 

1  o'  2  o  si  o  s2  o  1  o  2  o' 


‘2  “o' 

(37) 


the  corresponding  part  of  the  (continuous)  input  signal 
covariance  (33)  can  be  expressed  in  the  alternative  form 


"  %  exp(-|T|/TQ)  +  a*  exp(-|T|/T*)  = 


=  2  exp(-w^'-l)  [a^  cos(u^|t|)  +  sin(u)^lT|)] 


where  we  have  defined  the  real  and  imaginary  parts 


a  =  a  +  la-  , 
o  r  1  ' 


—  =  U)  +  10). 

^o 


(39) 


The  sampled  signal  covariance  for  component  (38)  is 


R°(kA)  =  2  exp(-a^|k|)  [a^  cos(a^|k|)  +  sin(a^|kl)]  ,  (40) 
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where  a^.  =  w^A  and  aj^  =  w^^A. 

The  spectral  component  associated  with  component  covariance 
(38)  is,  by  reference  to  (35)  and  (39), 


„  2  a  T_ 

Gg(f)  =  °  ° 


■C  + 


«  *  * 

2  a_  T 
o  o 


1  +  (w  Tq)^  1  +  (U)  t*)^ 


2  a  A  (a  +  ia.  )(a)  +  iu.  ) 

=  2  Re  -o — - — ^  =  4  Re  ^ 


2^1/2 
“  +  1/to 


+  (co^  +  iWj^)^ 


=  4 


a  (ji 
r  r 

f  2  2  21 

(*)^  +  (*)^  +  (0 

+  a .  (*)  ■ 

1  1 

(2^2  2'! 
+  (0^  -  (*) 

f  2  2) 

(0^  +  (w  -  w^)  J 

+  {Oi 

21 

+  (0.)^ 

.  “r  -  «i(‘"  -  “i)  .  .  “r  “r  “i(“  “i^ 

“  JL  o  o  ^  ^ 


wj  +  (w  -  w^)^ 


wj  +  (w  +  w^)^ 


(41) 


«i  i“r  «i  -  i“r  .  -  “i  i“r  .  ’  “i  '  ^“r 

X  -  '  \  +  '  - 1 -  X 


(0  +  0).  +10)  (0  +  {0.-  1(0^  (0-0).  +10) 
1  r  1  r  1 


0)  -  0)j^  -  10)^ 


If  this  spectral  component  Gg(f)  goes  negative  for  some  ranges  of 

f,  it  must  be  accompanied  by  additional  terms  in  summation  (35) 

to  keep  the  total  spectrum  G  (f)  nonnegative  for  all  frequencies 

s 

f.  spectral  component  (41)  has  four  poles  in  the  complex  o)-plane 

at  locations  co  =  (Oj^  ±  io)^  and  lo  =  -  o)^  ±  i(o^  . 

Even  if  coefficients  {a-J  and  A  • J  in  (35)  are  purely  real 

J  ®  J 

and  the  poles  are  distinct,  oscillatory  behavior  of  spectrum 
Gg(f)  is  still  attainable.  For  example,  with 

J  =  3,  (a^l  =  1,  8,  15,  (Oj)  =  11,  -48,  75,  (42) 
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the  spectrum  (35)  is  displayed  versus  wA  in  figure  6.  Actually 
plotted  is  the  scaled  spectrum 


^  °1  ‘1  _  , 

^  (UA)^  +  Aj 

where  we  have  used  (36)  to  eliminate  !fgjl*  negative 

coefficient  for  02  causes  the  dip  near  uA  »  4 . 


(43) 


Figure  6.  Spectrum  (43)  for  Parameters  (42) 
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INPUT  NOISE  CHARACTERIZATION 

The  presentation  in  this  subsection  exactly  parallels  that 
above  for  the  input  signal.  The  input  noise  covariance  is  again 
taken  to  be  a  sum  of  damped  exponentials  of  the  form 

M 

m=l 

where  the  M  (complex)  noise  time  constants  {r  are  all  distinct 

from  one  another.  (Some  of  them  may  equal  some  of  the  signal 

time  constants  {t  .)  in  (33),  if  desired.)  The  origin  value, 

s  J 

M  , 

R  (0)  =  P  =  o  ,  (45) 

m=l 

is  the  input  power  of  noise  sequence  n(kA)  in  (1). 

For  time  sampling  increment  A,  the  sampled  input  noise 
covariance  is,  from  (44), 

M 

R.,(kA)  =  exp(-b„lk|)  ;  b  =  —  .  (46) 

n  ‘ — r  m  m'  '  m  t  '  ' 

m=l  nm 

The  noise  exponential  parameters  are  distinct,  but  they  can 

be  complex;  however.  Re  b^^^  >  0.  (Some  of  the  (bj^^l  can  equal  some 

of  the  (aj)  in  (36).)  The  special  case  of  white  noise, 

2  2 
R^(kA)  =  6j^q,  is  handled  by  choosing  M  =  1,  -»  0, 

bj  +<». 
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FILTER  CHARACTERIZATION 


We  first  consider  an  analog  filter  with  a  rational  transfer 
function  with  N  distinct  poles,  namely 


1  +  i2nfT, 


(47) 


where  are  distinct  (complex)  time  constants.  H(f)  can  have 

up  to  N-1  zeros.  The  corresponding  impulse  response  is 


h(T  ) 


f  N  , 

HI  —  r" 

=  sn=l  ^hn  ^^hi 


for  T  i  0 


for  T  <  0, 


which  is  required  to  be  real  for  all  t.  The  sampled  impulse 
response  is  then  of  the  form 


A  h ( mA )  = 


2 _  tn  c  exp(-c  m)  for  m  ^  0  ,  c  =  - —  ,  (49) 

n=l  “  "  "  ^  ^hn 


where  the  N  filter  exponential  parameters  |Cj^)  are  distinct  from 
one  another;  also,  Re  c^  >  0. 

The  filter  correlation  function  is  available  by  substitution 
of  (49)  in  (10),  with  the  result 


'^n  “^n  ^n  ^  ' 


where 


N  c^ 

1  -  exp(- 


exp(-c^-c,^) 


for  1  1  n  1  N 


(51) 
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The  form  of  filter  correlation  ♦j^(jA)  in  (50)  is  identical  to  the 
sampled  input  signal  covariance  (36),  in  that  (50)  is  also 
composed  of  N  distinct  exponentials;  hence,  the  comments  and 
examples  given  there  are  relevant  here  also. 

A  useful  alternative  exists  to  the  calculation  of  constants 
via  (51).  Namely,  define  z  transform  filter 

00 

H(z)  =  YIZ  ^  h(mA)  .  (52) 

m=0 


Replace  n  by  k  in  (49)  and  substitute  the  result  into  (52); 
then  an  interchange  of  summations  yields 


H(z) 


N 


k=l 


^k 

1 

1  -  -  exp(-Cj^) 


(53) 


It  then  follows  immediately  by  comparison  with  (51)  that 

'’n  “  li(®*P(Cj^))  for  1  i  n  i  N  .  (54) 

When  H(z)  is  available  as  a  rational  function,  (54)  can  be  used 
directly,  instead  of  summation  (51). 

The  exponential  parameters  in  impulse  response  (49)  and 

corresponding  N-pole  filter  (53),  can  be  complex,  provided  that 
they  occur  in  complex  conjugate  pairs  and  that  their 
corresponding  coefficients  {+j^l  are  also  conjugate  pairs.  That 
is,  the  sum  of  two  conjugate  terms  of  the  form  of  (49)  must  be 
real  for  all  m  i  0.  This  generality  allows  for  oscillatory 
impulse  responses  such  as  yielded  by  narrowband  filters. 
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FILTER  OUTPUT  NOISE  COVARIANCE 

In  the  sequel,  we  will  find  use  for  the  following  closed  form 
for  the  discrete  convolution  of  two  exponentials; 


exp[-c|k|]  ®  exp[-b|k|]  s  ^  expt -c | j | -b | k- j | ]  = 

3 


sinh(c)  exp(-b|k|)  -  sinh(b)  exp(-c|k|) 


cosh(c)  -  cosh{b) 


[ixl  ^  ffsEif)] 


for  b  c 


for  b  =  c 


(55) 


The  filter  output  noise  covariance  is  given  by  (9).  If  we 
define 


C(kA)  =  R^(kA  -  jA)  ,  (56) 

then  the  stationarity  of  filter  output  noise  yj^(kA)  in  (8)  allows 
us  to  express  (9)  as 

Cj^(kjA,k2A)  *  C(k^A  -  k2A)  .  (57) 

Therefore,  we  can  concentrate  on  evaluation  of  C(kA)  in  (56)  for 
representative  filter  correlations  and  input  noise  covariances 

«n- 

In  particular,  we  will  use  the  general  filter  correlation 
(50)  and  input  noise  covariance  (46),  where  we  presume  that  none 
of  the  filter  parameters  (Cj^)  are  equal  to  any  of  the  input  noise 
parameters  That  is, 
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c  b  for  all  n,ni  . 
n  m 


(58) 


Substitution  of  (50)  and  (46)  in  (56)  and  use  of  (55)  results  in 


N 


M 


C(ka)  -  i;  JIZ  -l-n  '■n  exp(-c„|j|)  EZ  exp(-b„|k-j|)  = 

in=i 


:  n=l 


N  M 

>  !  >  \)„  6  72  exp[-c„|  j  |-b„|k-j  I  ]  = 

* — r  tz — r  ^n  n  n'^m^  n'-’’  m' 

n=  1  ni=  1  3 


M  N 


(59) 


m=l 


n=l 


where  we  defined  auxiliary  constants 


N  c  V  sinh(c  ) 

cosh(c^)  -  cosh(b^)  1  i  m  1  M  , 


M 


Y  =  U/  C  V 
'  n  T  n  n 


6  sinh(b_) 
'  m' 


n  ^n  n  n  r  cosh(b„)  -  cosh(c„) 
m=l  'in'  '  n' 


for  1  i  n  1  N 


(60) 


Condition  (58)  keeps  all  the  denominators  in  (60)  from  becoming 
zero.  On  the  other  hand,  if  one  of  the  filter  parameters  {c^}  is 
equal  to  one  of  the  noise  parameters  then  it  is  merely 

necessary  to  utilize  instead  the  second  line  of  (55)  for  that 
particular  n,m  pair  in  the  j  sum  in  the  middle  line  of  (59). 

End  result  (59)  for  the  filter  output  noise  covariance  is  a 
compact  expression  that  is  capable  of  being  quickly  computed, 
once  constants  {p^l  and  (y^)  have  been  evaluated  by  means  of 
(60)  and  stored.  Again,  we  encounter  a  sum  of  decaying 
exponentials,  albeit  of  M  +  N  terms  now. 
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FILTER  OUTPUT  NOISE  COVARIANCE  FOR  WHITE  NOISE  INPUT 

When  the  input  noise  is  white,  then  as  noted  under  (46),  we 
have  for  the  input  noise  covariance, 

Rn(kA)  -  ,  that  is,  M  -  1  ,  -»  +»  .  (61) 

2 

Use  of  this  result  in  (60)  yields  W,  0,  in 

'  *^1  'n  n^iinn 

which  case  the  filter  output  noise  covariance  (59)  becomes 

,  N 

C(kA)  ^  ^n  %  exp(-c^|k|)  ,  (62) 

n*l 

where  (v^)  are  given  by  (51). 

Actually,  (62)  is  a  special  case  of  the  general  white  noise 
result  obtained  by  substitution  of  (61)  into  (56);  namely,  for 
arbitrary  filter  A  h(mA),  the  filter  output  noise  covariance  is 

2  ... 

C(kA)  -  ^^(kA)  for  white  noise  input  ,  (63) 

where  filter  correlation  is  given  by  (10). 
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FILTER  OUTPUT  SIGNAL  COVARIANCE 

The  nonstationary  filter  output  signal  covariance  is  given  by 
double  sununation  (17)  in  conjunction  with  (18).  If  we  substitute 
input  signal  covariance  R  (kA)  given  by  (36),  along  with  filter 
impulse  response  A  h(mA)  given  by  (49),  into  (17),  there  follows 


In  order  to  evaluate  the  innermost  double  summation  in  (64), 
we  let  A  =  exp(aj),  B  =  exp(Cq),  C  =  exp(c^)  and  we  presume  that 
kj^  i  k2-  It  then  follows  from  (18)  that  Kj^  ^  ^2'  double 
sum  in  the  last  line  of  (64)  can  be  developed  thusly: 
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Ki  ^  K2 

=  JZZ  (C/A)"'  XZ  +  YU.  (CA)"*  ^  (B/A)P  -  (CB)"*  .  (66) 


Now  we  use  the  fact  that 


M  L  _  M+1 

ZZ  z  “  "  z  /  1  ,  L  ^  M  ,  (67) 

m=L  ^  ^ 


and  we  define  auxiliary  function 

k 

g(z,k)  2  2-2  z  1  .  (68) 

Then,  after  some  amount  of  manipulation,  the  sums  in  (66)  can  be 
expressed  in  the  closed  form 

2  2 
S(A,B,C,k3,Kj,K2)  =.  g(CB,k^)  g(CB,Kitl) 

+  g(B/A,K2+l)  [g(CA,K^+l)  -  g(CA,k^)]  -  g(BA,k^)  g(C/A,Kj+l)  , 

(69) 

provided  that  B  A  and  C  A.  (B  cannot  equal  1/A,  nor  can  C 
equal  1/A  or  1/B,  because  the  real  parts  of  c^  and  aj  are  always 
positive. ) 

We  now  employ  definition  (65)  in  order  to  express  the  filter 
output  signal  covariance  in  (64)  in  the  final  form 
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provided  that  none  of  the  filter  parameters  {c^^}  are  equal  to  any 
of  the  input  signal  parameters  {aj).  That  is, 

^  a^  for  all  n,j  .  (71) 

This  result,  (70),  cannot  be  further  reduced  without  the  various 
parameters  being  specified.  The  triple  summation  will  not  be 
overly  time  consuming  unless  the  input  signal  or  filter  are 
characterized  by  a  large  number  of  poles,  that  is,  large  J  or  N. 

The  result  in  (70)  holds  for  1  V.2)  the  range  for  >  V.2 
is  most  easily  handled  by  use  of  symmetry  relation  (15).  Also,  a 
significant  computational  aid  is  available  by  using  the 
recurrence  relation  for  the  g(z,k)  function  in  (68),  namely 
g(z,k)  =  z  g(z,k-l) . 

In  summary,  k^,  kj^,  k^,  k^  are  given  integers  satisfying 

^  k.  ,  k^  1  k. ,  k,  i  k. .  Also,  =  max(k^,k^), 
a  b  e  d  a  d  o  'a  c 

Kj  =  min( kj^ , kj^) ,  K2  =  min(k2,kj^).  Since  we  keep  k^  <;  k2,  then 
Kj  i  K2.  Integers  kj^  and  k2  must  vary  over  the  range 

1  k^  i  k2  i  k^  ,  (72) 

meaning  that  varies  in  the  range  min(K^,kj^)  to  min(k^,kjj). 

From  (69),  we  see  that  we  have  to  compute  function  values 
g(CB,k^),  g(CA,k^),  g(BA,k^),  as  well  as  the  arrays  of  function 
values  of  g(CB,k),  g(CA,k),  g(B/A,k),  g(C/A,k)  for  k  in  the  range 
min(K^,kj^)  +  l  to  min(  k^,  kj^ ) +  1 . 
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PROGRAMS 

A  program  listing  for  the  white  noise  input  case  is  presented 
in  appendix  D;  it  consists  of  two  separate  parts.  The  first  part 
computes  the  filter  output  signal  and  noise  covariance  matrices 
for  unity  input  signal  and  noise  powers,  the  total  filter  output 
covariance  matrix  of  elements  C^,  in  (23)  for  the  various  signal- 
to-noise  ratios  of  interest,  the  weighted  covariance  matrix  C  of 
elements  C,  in  (26)  and  its  corresponding  eigenvalue  matrix  A  in 

a 

(27).  The  output  of  these  eigenvalues  to  storage  completes  the 
first  part. 

The  second  part  lists  the  program  that  takes  in  these 
eigenvalues  and  computes  the  exceedance  distribution  functions, 
for  noise-alone  as  well  as  for  signal  present,  from 
characteristic  function  (32);  these  are  the  false  alarm  and 
detection  probabilities,  respectively.  The  precise  method  of 
handling  all  the  widely  different  signal-to-noise  ratios  of 
interest  is  presented  in  the  next  section.  At  the  end  of 
appendix  D,  there  is  also  a  listing  of  the  simulation  program 
utilized  to  verify  these  theoretical  derivations  and  results. 

The  corresponding  program  for  the  colored  noise  case  is 
listed  in  appendix  E.  A  similar  breakdown  into  two  parts  has 
been  adopted;  however,  once  the  eigenvalues  have  been  computed  at 
the  end  of  the  first  part,  the  identical  program  in  part  2  of 
appendix  D  is  used  for  distribution  calculations  and  is  therefore 
not  listed  again.  Also,  the  simulation  program  use  to  check  the 
colored  noise  case  is  similar  to  that  in  part  3  of  appendix  D. 
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ON  CALCULATION  OF  RECEIVER  OPERATING  CHARACTERISTICS 

The  evaluation  of  the  receiver  operating  characteristics  for 
a  number  of  different  input  signal-to-noise  ratios  typically 
involves  the  Fourier  transforms  of  a  set  of  characteristic 
functions  of  the  decision  variable,  which  frequently  have  widely 
different  ranges  and  rates  of  variation.  In  order  to  make  these 
calculations  tractable  and  easily  plotted,  it  is  necessary  to 
resort  to  FFTs  [2]  with  carefully  chosen  sizes  and  sampling 
increments.  We  now  present  the  reasoning  behind  our  choices  of 
these  parameters  and  their  interrelationships. 

For  ease  of  presentation,  we  will  consider  the  numerical 
evaluation  of  a  probability  density  function  Pq(u)  from  a  given 
(noise-only)  characteristic  function  the  extension  to  an 

exceedance  distribution  function  [2]  is  immediate.  We  have 

Po^^^  "  Zn  I  exp(-if;u)  f^(^)  = 

=  A  Y]  exp(-ikA^u)  f  (kA  )  s  p  (u)  ,  (73) 

2n  o  ^  o  o  o  ^o 

where  the  trapezoidal  approximation  with  sampling  increment  A^ 
in  I  has  been  employed.  The  approximation  Pq(u)  defined  by  (73) 
is  an  aliased  version  of  Pq{u)  and  has  period  u^  =  2n/L^  in  u. 

In  order  to  keep  aliasing  effects  negligible  in  Pq(u),  it  will  be 
necessary  to  choose  I  increment  A^  small  enough  that  the  aliasing 
lobes  at  separation  u^  don't  overlap.  Therefore,  we  can  restrict 
the  evaluation  of  Pq(u)  to  N^  equally  spaced  points  over  that 
interval  u^,  according  to 
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Po(i^)  =  2i  ^  exp(-i2nkn/N^)  for  0  1  n  1  N^-1  . 

O  O  K 

(74) 

Furthermore,  by  collapsing  the  sequence  {f^(kA^))  modulo  N^, 

(74)  can  be  accomplished  exactly  and  efficiently  as  an  N^-point 
FFT  when  Nq  is  a  power  of  2;  see  [2].  The  increment  in  argument 
u  of  Pq(u)  in  (74)  and  the  maximum  useful  value  of  u  are 


& 


o 


2n 

N  A  ' 
o  o 


(75) 


Now  suppose  that  we  also  want  to  evaluate  the  probability 
density  function  P2^(u)  corresponding  to  a  different  (signal 
present)  characteristic  function  according  to 

Pl(u)  =  ^  J  d?;  exp(-i?;u)  fj(j;)  = 

=  2n  E  exp(-ikAjU)  fj(kAj)  =  pj(u)  ,  (76) 

k 

where  the  sampling  increment  in  I  is  now  A^.  By  identical 
reasoning  to  that  above,  we  can  get  samples  of  periodic 
(aliased)  approximation  Pj(u)  according  to 

Pi(f^)  “  2n  ^1  E!  exp(-i2nkn/Nj)  fj(kAj)  for  0  i  n  i  N^-1  , 
11  k 

(77) 

which  can  be  realized  as  an  N^-point  FFT.  The  increment  in  u  and 
the  maximum  useful  value  of  u  for  this  latter  case  are 


& 


1 


2n 

NiAj 


(78) 
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If  we  now  want  to  eliminate  the  arguments  of  functions  Pq(u) 
and  P2(u)  and  be  able  to  directly  plot  Pj^Cu)  versus  Pq{u)  (that 
is,  y  versus  x),  this  can  be  easily  accomplished  if  we  force  the 
u  increments  in  (74)  and  (77)  to  be  equal,  namely,  =  6^.  But, 
then,  according  to  (75)  and  (78),  this  means  choosing 

"l''!  *  Vo  • 

Since  integers  and  will  be  limited  to  powers  of  2,  this 
will  mean  repeated  halving  of  the  K  increments  and  Aj^  in  (73) 
and  (76),  respectively.  The  values  of  Pj^  yielded  by  (77)  can 
then  be  directly  plotted  versus  those  values  of  p^  yielded  by 
(74),  up  to  n  =  min(NQ,Nj)  =  N^. 

There  are  several  alternative  procedures  that  could  have  been 
adopted.  For  example,  one  option  is  to  halve  the  K  increment, 
that  is,  take  Aj^  =  keep  Nj^  =  N^.  This  will  require 

discarding  every  other  value  put  out  by  FFT  (74),  in  order  that 
the  remaining  values  occur  at  the  same  u  arguments  yielded  by 
(77).  This  is  wasteful  of  FFT  (74)  and  has  not  been  adopted 
here.  Nor  have  we  employed  the  possibility  of  interpolation  of 
(77)  to  determine  approximately  what  the  values  of  Pj^  would  be  at 
the  arguments  of  in  (74). 

The  need  to  decrease  the  size  of  the  J;  increment  from  A^  in 
(74)  to  Aj  in  (77)  is  fundamental.  It  arises  from  the  fact  that 
as  the  input  signal-to-noise  ratio  is  increased,  the  contribution 
of  the  corresponding  probability  density  function  Pj^(u), 
k  =  0,1,2,...,  moves  to  higher  values  on  the  u  scale  (thereby 


33 


TR  8913 


leading  to  the  desired  higher  detection  probabilities).  In  order 
that  approximation  not  be  severely  aliased,  the  K  increment 

must  therefore  be  decreased;  for  example,  see  upper  limits  u^ 
and  Uj^  in  (75)  and  (78),  respectively. 

With  these  points  in  mind,  the  following  procedure  has  been 
adopted  and  utilized  in  the  programs  in  this  report.  First,  for 
the  noise-only  probability  density  function  p^{u) ,  a  satisfactory 
value  for  i;-increment  in  (73)  is  found  such  that  aliasing  is 
insignificant  in  Pq(u).  This  selection  of  A^  is  arrived  at  by 
looking  at  a  plot  of  (74)  and  modifying  A^  appropriately  by  trial 
and  error.  An  FFT  size  of  *  128  is  utilized  to  evaluate  (74), 
which  is  then  stored;  this  size  for  has  generally  proven  to  be 
sufficient  to  keep  track  of  the  variations  of  p^(u) , 

When  we  encounter  the  characteristic  function  fj^(JI)  and 
probability  density  function  Pj^(u)  for  the  first  (lowest)  nonzero 
signal-to-noise  ratio  of  interest,  it  is  usually  necessary  to 
decrease  the  K  increment  to  Aj^  »  order  to  control 

the  aliasing  inherent  in  approximation  pj^(u).  At  the  same  time, 
is  scaled  up  by  2  or  4,  according  to  (79),  thereby  maintaining 
the  same  u  arguments  for  (74)  and  (77),  as  desired.  Again,  a 
plot  of  aliased  version  pj(u),  obtained  by  means  of  (77),  serves 
as  the  decision  point  for  making  an  acceptable  choice  of  value 
for  Aj^. 

For  the  successively  larger  signal-to-noise  ratios  and  their 
corresponding  probability  density  functions  Pj^(u),  k“l,2,..., 
occasional  repeated  halvings  of  the  £;  increment  according  to 
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Aj^  =  made,  but  only  when  needed  to  keep  aliasing 

insignificant  in  pj^(u).  At  those  times  that  such  a  halving 
is  necessary  and  made,  the  FFT  size  is  doubled  according  to 
Nr  =  2  Otherwise,  Aj^  and  Nj^  are  kept  at  their  same 

values  as  used  for  the  k-1  run.  By  the  time  the  largest 
signal-to-noise  ratio  case  of  interest  is  reached,  the  FFT  size 
Nj^  can  get  rather  large.  We  are  limited  in  our  application  to  a 
maximum  FFT  size  of  16384;  if  aliasing  is  significant  at  this 
size,  we  are  unable  to  accurately  numerically  handle  these  larger 
signal-to-noise  ratios  by  the  adopted  procedure.  Interpolation 
would  then  be  the  recommended  alternative,  as  alluded  to  above. 

The  procedure  above  has  been  explained  in  terms  of 
probability  density  functions  rather  than  the  exceedance 
distribution  functions  which  are  of  actual  interest.  That  is, 
the  relevant  exceedance  distribution  functions  are  plots  of  the 
detection  probabilities  for  various  signal-to-noise  ratios  versus 
the  false  alarm  probability.  The  only  basic  change  is  to  replace 
the  top  lines  of  (73)  and  (76)  by  the  exceedance  distributions, 
which  are  also  obtained  by  means  of  Fourier  transforms,  namely 
[2;  (5)  -  (6)] 


Ek(u) 


+00 


0+ 


exp(-iJ;u) 


fkiH 


for  k  *  0,1,2,... 


(80) 
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RESULTS 


COMPARISON  WITH  SIMULATION 

In  figure  7,  the  cumulative  (CDF)  and  exceedance  distribution 
functions  (EDF)  for  a  white  noise  input  excitation  are  presented. 
The  input  signal  and  filter  are  both  characterized  by  one  pole, 
that  is,  J  =  1  in  (33)  -  (36)  and  N  »  1  in  (47)  -  (50).  The 
simulation  employed  1E6  trials  and  verifies  the  theoretical 
calculations  down  to  the  probability  level  of  lE-5  plotted.  A 
listing  of  the  simulation  program,  including  all  the  parameter 
values  that  were  utilized  here,  is  given  at  the  end  of  appendix 
D.  The  gating  and  observation  parameters  are  =  4,  =  11  and 

k^  *  2,  k^  =  16,  respectively;  thus,  it  can  be  seen  that  the 
transient  buildup  and  decay  are  a  prominent  part  of  the  filter 
output  signal  for  this  particular  example. 

The  results  in  figure  8  pertain  to  a  colored  noise  input  with 
one  pole,  that  is,  M  =  1  in  (44)  -  (46).  Again,  1E6  trials  were 
used  in  the  simulation  and  they  confirm  the  theoretical  result 
down  to  the  lE-5  level  of  probabilities.  A  listing  of  the 
simulation  program  is  given  at  the  end  of  appendix  D. 

Absolute  time  is  unimportant  to  the  performance  of  the  energy 
detector  in  figure  1.  That  is,  any  common  constant  can  be  added 
to  or  subtracted  from  k^,  kj^,  k^,  k^,  without  changing  the 
receiver  operating  characteristics.  Thus,  k^  could  always  be 
selected  as  0  without  loss  of  generality. 
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Figure  7 .  Cumulative  and  Exceedance  Distributions  for  White  Noise 


Figure  8.  Cumulative  &  Exceedance  Distributions  for  Colored  Noise 
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RECEIVER  OPERATING  CHARACTERISTICS 

The  receiver  operating  characteristics,  that  is,  detecticr 
probability  versus  false  alarm  probability  Pp  for  a  set  of 
signal-to-noise  ratios,  for  the  system  in  figure  1  will  be 
considered  in  this  section.  The  first  example  is  evaluated  for 
the  following  set  of  parameters,  where  all  times  are  in  seconds. 
Also,  all  parameters  are  real. 

k^  =  4,  11»  k^  =  2,  k^»16,  A  =  .2, 


J  =  3, 

(ajl 

=  {11 

-48 

75),  {T^jl  = 

{1  1/8 

1/15), 

M  =  2, 

IPml 

=  {39 

60), 

t^nml  = 

.4), 

11 

{^'nl 

=  {1 

2  3 

l^hnJ  =  1 

.3  .5 

.7  .9). 

(81) 

Twenty  nonzero  values  of  the  system  input  signal-to-noise  ratio 

2 

R  -  ^  (82) 

°n 

were  utilized,  namely  R  =  5(1.2)27.8  dB.  For  computation  of  the 
exceedance  distribution  function  Pp  directly  from  the 
characteristic  function,  the  initial  I,  increment  in  (73)  and 
(74)  was  taken  as  .00025;  this  yielded  round-off  errors  less  than 
lE-15  in  the  false  alarm  probability  calculation.  Repeated 
occasional  halving  of  the  t  increment,  as  explained  in  the 
previous  section,  was  utilized,  eventually  requiring  an  FFT  size 
of  16384  for  the  largest  signal-to-noise  ratios  cases  for  P^. 

The  receiver  operating  characteristics  are  plotted  on  normal 
probability  paper  in  figure  9  and  cover  a  wide  range  of  values. 
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ranging  from  lE-10  for  Pp  to  .999  for  P^.  Since  Gaussian  random 
variables  would  plot  on  such  paper  as  a  set  of  parallel  straighL 
lines,  the  curvatures  of  these  results  illustrate  that  the 
Gaussian  assumption  for  decision  variable  z  in  figure  1  is  not 
warranted,  at  least  for  this  example.  The  major  reason  for  this 
behavior  is  the  relatively  small  number  of  samples  used  in  the 
detector;  see  the  top  line  of  (81).  This  is  also  a  partial 
explanation  for  the  rather  large  values  of  the  per-sample 
signal-to-noise  ratio  R  required  at  the  high  quality  region  at 
the  top  left  of  figure  9.  Another  factor  to  notice  is  that  the 
input  signal  duration  is  only  kj^A  -  k^A  =  1.4  seconds,  meaning 
that  the  filter  output  signal  never  reaches  steady  state  before 
the  input  signal  is  turned  off;  all  these  transient  signal 
effects  and  their  limitations  on  performance  have  been  included 
exactly  in  the  analysis  and  numerical  results  presented  here. 

The  effect  of  halving  the  sampling  increment  A  in  figure  1  is 
considered  in  the  next  example  but  done  in  such  a  manner  as  to 
keep  the  total  input  signal  duration  the  same.  That  is,  the 
parameter  values  in  (81)  are  kept  the  same  except  for  the 
following  changes: 

ka  =  8,  kjj  =  22,  k^  =  4,  k^  =  32,  A  =  .1.  (83) 

Notice  that  all  absolute  times,  such  as  k_A,  are  kept  fixed.  The 

O 

corresponding  receiver  operating  characteristics  are  plotted  in 
figure  10.  Performance  has  been  degraded  by  a  couple  of  dB.  For 
example,  to  realize  Pp  =  lE-3  and  P^  =  .5,  the  per-sample  input 
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signal-to-noise  ratio  R  must  now  be  16  dB,  whereas  it  was  13.8  dB 
in  the  previous  figure.  This  is  a  degradation  of  2.2  dB.  The 
main  reason  for  this  behavior  is  again  the  inability  of  the 
filter  output  signal  to  reach  steady  state  and  thereby  contribute 
significantly  to  the  summer  output.  By  contrast,  the  filter 
output  noise  is  in  steady  state  (by  assumption)  for  both 


examples . 

We  now  return  to  the  original  sampling  increment  A  =  .2 
seconds  employed  in  (81).  When  the  observation  interval 
coincides  with  the  input  signal  nonzero-excitation  duration,  that 
is,  =  4,  ~  ^d  ~  receiver  operating 

characteristics  in  figure  11  illustrate  additional  improvement. 
For  example,  probabilities  Pj,  =  lE-3  and  Pj^  =  .5  can  now  be 
realized  with  R  *  13.1  dB,  which  is  an  improvement  of  .7  dB 
relative  to  the  signal-to-noise  ratio  required  for  the  broader 
observation  interval,  k^  =  2,  k^  =  16,  used  in  (81). 

Finally,  an  example  with  a  wide  observation  interval,  namely 
k^  =  0,  k^  =  25,  is  displayed  in  figure  12.  The  desired 
probabilities  can  now  be  achieved  only  if  the  input  signal-to- 
noise  ratio  is  increased  to  14.7  dB,  a  degradation  of  1.6  dB 
relative  to  the  best  case  above. 


These  examples  illustrate  the  utility  of  being  able  to 
investigate  quantitatively  and  accurately  the  effects  of 
nonstationarity  in  a  system,  without  having  to  make  questionable 
assumptions  about,  for  example,  how  closely  steady  state  was  or 
was  not  realized.  They  also  afford  a  dependable  verification  or 
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rejection  of  the  Gaussian  approximation  for  the  decision 
variable;  in  a  related  work  [3]/  the  latter  approximation  was 
found  CO  be  too  optimistic  ror  most  working  ranges  of  this 
detection  system. 
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Figure  10.  Receiver  Operating  Characteristics;  Example  B 
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SUMMARY 

Programs  have  been  written  which  enable  exact  analysis  of 
performance  of  the  mismatched  bandpass  energy  detector  of  figure 
1,  both  for  the  white  noise  input  case  as  well  as  the  colored 
noise  input  case.  These  results  allow  for  arbitrary  sampling 
time  increment  A  and  for  arbitrary  input  signal  spectra,  input 
noise  spectra,  and  filter  transfer  functions.  By  these  means, 
exact  quantitative  evciluations  and  degradations  can  be  determined 
for  various  combinations  of  uncertainty  regarding  the  input 
signal  time  location  and  duration  as  well  as  its  center 
frequency,  bandwidth,  and  spectrum.  Included  in  the  analysis  and 
programs  are  the  buildup  and  decay  (or  any  portion  thereof)  of 
the  nonstationary  output  signal  from  the  filter,  when  excited  by 
a  burst-like  input  signal,  regardless  of  the  lengths  or  locations 
of  the  excitation  and  observation  intervals.  Both  programs  have 
been  compared  with  simulation  results  and  confirmed  down  to  as 
low  a  probability  level  as  possible  consistent  with  the  number  of 
trials  used. 
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APPENDIX  A.  NUMERICAL  EVALUATION  OF  NOISE  COVARIANCE  R^^ 

If  the  noise  covariance  Rjj(t)  is  not  available  in  analytic 
form,  but  the  noise  spectrum  G^(f)  is  specified,  the  following 
numerical  procedure  can  be  employed  to  evaluate  the  required 
samples  R^(kA).  We  begin  with 

Rn(T)  =  J  df  exp(i2nfT)  G^(f)  = 


S  Aj  XI  exp(i2nmA^T)  G^(mAj)  s  '  (A-1) 

m 

where  the  trapezoidal  rule  with  frequency  increment  was  used. 
But  in  (A-1)  can  be  developed  as 


^n^^^  “  J  exp(i2nfT)  G^(f)  A^ 

=  R„(r)  ®  -  E  R„(t  -  • 


(A-2) 


Now  suppose  that  the  noise  covariance 


R^(t)  =0  for  |t|  >  Tq  . 


(A-3) 


Then,  in  order  to  avoid  significant  aliasing  in  (A-2),  we  must 
take  A^  small  enough  that 

1 


257  >  '^o  ' 


(A-4) 


in  which  case  (A-2)  yields 


Rn(T)  =  R„(r)  for  |t|  <  ^ 


(A-5) 
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The  required  samples  of  follow  from  (A-l)  according  to 

R^(kA)  =  Aj  X]  exp{i2nmAjkA)  G^(mAj)  .  (A-6) 

m 


Now  take  frequency  increment 


(A-7) 


where,  according  to  (A-4),  integer  N  must  be  taken  large  enough 
that  NA/2  >  Tq,  that  is, 


N  > 


(A-8) 


Then,  from  (A-5)  -  (A-7), 

Rn(^^)  “  *  nI  ^  exp(i2nmk/N)  for  |k|  <  |  .  (A-9) 

m  '  ' 


Alternatively,  this  can  be  expressed  as 


N-1 


R^(kA)  £  Rj^(kA)  =  IZ  exp(i2nmk/N)  for  |k|  <  |  ,(A-10) 

m=0  ' 


where  collapsed  (prealiased)  noise  spectiu^a  and  sequence 


G„(() 


"  C  G„(f  -  i)  ,  .  J:  oJjjE  -  1)  tor  0  i  m  S  N-l  . 

(A-ll) 


The  procedure  in  (A-10)  can  be  efficiently  realized  as  an  N-point 
FFT  of  N  nonzero  numbers.  The  major  conditions  that  must  be  met 
are  (A-8)  in  conjunction  with  (A-3). 
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APPENDIX  B.  NUMERICAL  EVALUATION  OF  FILTER  CORRELATION 
Suppose  that  we  want  digital  filter  correlation 

E  h(inA)  h(mA  -  jA)  ,  (B-1) 

m 

defined  in  (10),  but  all  that  we  have  is  the  filter  transfer 
function  H(f),  where 

h(T)  =  J  df  exp(i2nfT)  H(f)  for  all  x  .  (B-2) 

Impulse  response  h(T)  is  presumed  real.  The  direct  relationship 
between  ♦^^{jA)  and  H(f)  is  the  subject  of  this  appendix. 

First,  define  periodic  function 

H(f)  =  E  H(f  -  f)  for  all  f  .  (B-3) 

m  '  ^ 

Then,  there  follows 


H(f)  =  H(f)  e  » 


=  J  dr  exp(-i2nfT)  h(T)  A 


-  E 

m 


exp(-i2nfAm)  A  h(mA)  ,  (B-4) 


where  we  used  the  fact  that  convolution  in  the  frequency  domain 
is  equivalent  to  Fourier  transformation  of  a  product  in  the  time 
domain.  That  is,  H(f)  is  the  Fourier  transform  of  the  samples  of 
impulse  response  h(T)  taken  at  time  increment  A.  However,  care 
must  be  taken  at  points  of  discontinuity  of  h(T);  for  example,  if 
h(T)  is  discontinuous  at  kA,  the  contribution  to  the  right-hand 
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side  of  (B-4),  at  m  =  k,  is  i[h(kA+)  +  h(kA-)].  This  point  and 
an  example  are  discussed  more  fully  in  appendix  C. 

Now,  there  follows  from  (B-1)  and  (B-2), 

♦h(j^)  *  H  J  df  exp(i2nfAm)  H(f)  J  du  exp[-i2nuA(m-j ) ]  H*(u)» 
m 

=  A  JJ  df  du  H(f)  H*(u)  exp(i2nuAj)  A  Y1  exp[i2n(f  -  u)nm]  * 

m 

=  A  JJ  df  du  H(f)  H*(u)  exp(i2nuAj)  “ 

=  A  ^  J  df  H(f)  H*(f  -  I)  exp[i2n(f  -  ^)Aj]  - 

=  A  J  df  exp(i2iifAj)  H(f)  XI  ■  ?)  • 

''  m  '  ' 


=  A  J  df  exp(i2nfAj)  H(f)  H*(f)  , 


(B-5) 


where  we  used  (B-3). 

Now,  let  denote  the  following  interval  of  length  1/A  on 
the  f  axis; 


(B-6) 


Then,  the  filter  correlation  follows  from  (B-5)  according  to 

=  A  XI  J  df  exp(i2nfAj)  H(f)  H*(f)  » 
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-A  E  J  du  exp[i2n[u  +  fj^j]  **(“■•■?)  «*(“)  = 

=  A  J  du  exp(i2nuAj)  H*(u)  ]E  '*’?)“ 

.5/A 

-A  [  du  exp(i2nuAj)  |H(u)|^  ,  {B-7) 

-.5/A 


where  we  let  u  *  f  -  n/A  and  used  (B-3)  and  the  periodicity  of 
H(f).  This  result  indicates  that  we  must  first  alias  the  given 
transfer  function  H(f)  according  to  (B-3)  and  then  Fourier 
transform  its  magnitude-square  over  an  interval  of  length  1/A. 

An  alternative  approach  that  leads  to  the  same  result  (B-7) 
is  to  use  (B-2)  and  (B-3)  in  the  form 

h(mA)  *  J  df  exp(i2iifAm)  H(f)  *  J  df  exp(i2nfAm)  H(f)  .  (B-8) 

^o 

Use  of  this  latter  result  in  (B-1)  leads  to  filter  correlation 


♦jj(j^)  “  E  J  df  exp(i2nfAm)  H(f)  J  du  exp[-i2nuA(m-j )  ]  H*(u)  = 
m  , 


du  exp(i2nuAj)  H(f)  H*(u)  A  Y1  exp[i2n(f  -  u)Am] 

m 

df  du  exp(i2nuAj)  H(f)  H*(u)  E!  “  z]  * 
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=  A  J  df  exp{i2ii£Aj)  lH{f)|^  ,  (B-9) 

since  only  one  impulse,  at  u  -  f,  lies  in  interval  1^.  Result 
(B-9)  is  identical  to  (B-7).  Notice  that  we  do  not  get 


A  J  df  exp{i2nfAj)  |H(f)|^  »  A  J  df  exp{i2nfAj)  Y.  ”[^'^5)1 


(B-10) 


As  for  the  actual  numerical  evaluation  of  {B-9),  suppose  we 
sample  the  integral  on  f  with  increment  A^  *  1/{KA)  and  use  the 
trapezoidal  approximation;  then  there  follows 


♦ 


K/2 


k=-K/2 


exp 


(B-U) 


where  weights  {Wj^}  are  given  by 


for  !k|  -  K/2 
U  for  |k|  <  K/2 


(B-12) 


By  using  the  periodicity  K  of  the  exp  and  H  terms  in  (B-11),  that 
sum  may  be  written  exactly  in  the  standard  FFT  form 


♦h(3A) 


I 

K 


K-1 

>  '  exp(  i2ix  jk/K) 

k=0 


for  |j|  <  I  . 


(B-13) 


This  is  a  good  approximation,  provided  that  the  FFT  size  K 
satisfies  K  >  2tj^/A,  where  is  the  effective  duration  of  filter 
correlation 


54 


TR  8913 


APPENDIX  C.  DISPLACED  SAMPLING  AND 
FOURIER  TRANSFORM  OF  DISCONTINUOUS  FUNCTION 

GENERAL  RELATIONS 

Let  a(t)  and  A(f)  be  a  Fourier  transform  pair: 

A{f)  =  J  dt  exp(-i2nft)  a(t)  .  (C-1) 

Also,  let  b(t)  and  B(f)  be  a  Fourier  transform  pair.  Then, 
Parseval's  theorem  states  that 

J  dt  a(t)  b*(t)  “  J  df  A(f)  B*(f)  .  (C-2) 

We  now  apply  this  relation  to  the  case  where 

b(t)  *  A  6^(t  -  a)  ,  B(f)  -  exp(-i2nfa)  .  (C-3) 

There  follows 

L  s  A  X]  a{nA  +  a)  »  f  dt  a(t)  A  5.(t  -  a)  *  (C-4) 

n  ^ 

=  J  df  A{f)  exp(i2nfa)  =  Y1  exp{i2nna/A)  s  r  .  (C-5) 

If  the  samples  required  in  the  sums  in  {C-4)  and  (C-5)  encounter 
a  discontinuity  of  a(t)  or  A(f),  the  value  used  in  (C-4)  or  (C-5) 
must  be  taken  as  the  average  value  approached  from  both  sides  of 
the  discontinuity.  Also,  the  sums  in  (C-4)  and  (C-5)  may  have  to 
be  interpreted  as  principal  value,  if  necessary  for  convergence. 
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EXAMPLE 


a(t) 


^exp[-(b+ic)t] 
[  0 


for  t  >  0 
for  t  <  0 


b  >  0  . 


(C-6) 


A(f) 


1 _ 

b  +  i(c  +  2nf) 


for  all  f  . 


(C-7) 


Function  a(t)  is  discontinuous  at  t  =  0,  while  A(f)  is  continuous 
for  all  f. 

Let  0  <  a  <  A;  then  the  sum  in  (C-4)  is 


L  *  A  >  '  exp[-(b+ic) (nA+a) 1 

n=0 


A  expr -a(b+ic)  1 
1  -  exp[ -A{b+ic ) ] 


for  0  <  a  <  A  . 

(C-8) 


At  the  same  'ime,  the  sura  in  {C-5>  is 


R 


E 


exp( i2nna/A ) 
b  +  i{c  +  2nn/A) 


for  all  a 


{C-9) 


When  a  /  0  (or  an  integer  multiple  of  A),  the  phase  factor  in  the 
numerator  of  (C-9)  yields  a  convergent  sum  for  the  real  and 
imaginary  parts,  without  the  need  for  a  principal  value 
interpretation.  It  has  also  been  verified  numerically  that  (C-8) 
and  (C-9)  are  equal,  as  predicted  by  (C-4)  and  (C-5);  that  is. 


Eexp( i2nna/A) 
b  +  i(c  +  2nn/A) 

T1  '  ' 


A  expf -a(b->-ic )  1 
1  -  exp[-A(b+ic) ] 


for  0  <  a  <  A 


(C-10) 
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On  the  other  hand,  if  a  *  0,  then  the  sum  in  (C-4)  is, 
accounting  for  the  discontinuity  of  a(t). 


CO 

I  +  A  YZZ  exp[-(b+ic)nA]  »  ^  y~ 
n*l 


A  1  +  expr-A(b-t-ic)  1 
2  1  -  exp[-A(b+ic) ] 


(C-11) 


Also,  then,  the  sum  in  (C-5)  is 


_ 1 

“  b  +  i(c  + 


b  -  i(c  +  2nn/A] 


(c  +  2nn/A)  n  b^  +  (c  +  2iin/A)^ 


(C-12) 


the  imaginary  part  of  which  must  be  interpreted  as  a  principal 
value  sum.  Again,  it  has  been  numerically  verified  (next  page) 
that  (C-11)  and  (C-12)  are  equal;  that  is. 


y  _ , _ L 

4-  b  +  i(c  -T 


1 _  ^  A  1  -f  expr-A(b+ic)  1 

(c  T  2nn/A)  *  2  1  -  exp[-A(b+xc) ] 


(C-13) 


Notice  that  the  limit  of  (C-10)  as  a  0+  does  not  yield  (C-13). 

If  we  apply  result  (C-11)  to  H(f)  in  the  last  entry  in  (B-4), 
we  have,  for  example  (304), 


5/<:v  _  a  1  +  exp(-a-i2nfA) 

2  1  -  exp(-a-i2nfA)  ’ 


(C-14) 


By  contrast,  the  limit  of  (C-10)  as  a  ->  0+  yields 


1  -  exp( -a-i2nf A ) 


(C-15) 
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BASIC  PROGRAM  FOR  NUMERICAL  VERIFICATION 


Bs=.57  !  b 

Cs=.71  I  c 

De=.93  !  delta 

fll=.l  i  alpha 

PRINT  6s,Cs,Ile,R1 
CfiLL  ExpC-fll«Bs,-fll*Cs,Nr,Ni ) 

CALL  Exp< -De*Bs , -De*Cs , Er , E i ) 

CALL  Div<De*Nr,De»Ni , 1-Er,-Ei ,Lr,Li ) 
IF  A1=0.  THEN  Lr*=Lr-.5*De 


PRINT  Lr.Li 
A=2*PI/De 
B=R*A1 
B2=Bs*Bs 
Rr=Ri =0. 

DOUBLE  Ns 

FOR  Hs=-1E5  TO  1E5 

T=B*Ns 

C=COS<T) 

S=SIN<T) 

T=Cs+R*Ns 
D=i:2iTr. 
Dr=<C*Bs+S*T)^D 
Di 9<S*Bs-C*T>/D 
Rr=Rr+Dr 
Ri =Ri +Di 
NEXT  Ns 
PRINT  Rr,Ri 
END 


LEFT-HAND  SIDE  L 


INTEGER 

PRINCIPAL  VALUE  SUM 


RIGHT-HAND  SIDE  R 


SUB  DiwCXl, Yl,X2, Y2,U,V) 

T=X2*X2+Y2*Y2 

U=<X1*X2+Y1*Y2>/T 

V=<Y1*X2-X1*Y2)/'T 

SUBEND 


Z1/'Z2 


SUB  Exp<X,Y,U,V) 
E=EXP<X) 
U=E*COS<Y) 
V=E*SIN< Y) 

SUBEND 


EXP<Z> 


.  57 

. 729357897734 
. 729357647983 


.71 

-.80564044434 

-.805640755434 


.  57 

1 . 19310532814 
1 . 1935371 1305 


.71 

.806028668009 

-.806028668087 


.  57 

1 . 07135515262 
1 . 07135537888 


.71 

.839119622089 

-.839119622093 
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APPENDIX  D.  PROGRAMS  FOR  WHITE  INPUT  NOISE 

Three  BASIC  programs  are  listed  in  this  appendix.  All  the 
signal,  noise,  and  filter  parameters  are  restricted  to  be  real 
here.  The  first  program  computes  the  filter  output  signal  and 
noise  covariance  matrices  for  a  white  noise  input,  the  sums  of 
these  matrices  scaled  by  the  various  signal-to-noise  ratios  and 
then  weighted,  and  the  eigenvalues  which  are  then  stored.  The 
particular  subroutine  listed  here,  SUB  Svd,  actually  computes  the 
eigenvectors  in  addition;  a  faster  alternative  would  be  to 
replace  this  subroutine  by  one  which  calculates  only  the 
eigenvalues.  Twenty  different  nonzero  values  of  the  input 
signal-to-noise  ratio  are  allowed  in  the  program  and  must  be 
chosen  and  input  by  the  user. 

These  eigenvalues  serve  as  the  input  to  the  second  program 
which  computes  the  cumulative  and  exceedance  distribution 
functions  according  to  the  method  given  in  [2]  and  then  plots 
the  receiver  operating  characteristics  by  elimination  of  the 
threshold  variable  according  to  the  method  described  in 
(73)  -  (80)  in  the  main  text. 

The  third  program  is  the  simulation  program  used  to  check  the 
two  programs  above.  It  is  written  for  single-pole  processes  and 
allows  for  a  colored  input  noise  process  but  could  be  easily 
modified  to  handle  more  general  processes.  A  run  of  1E6  trials 
took  7  hours  on  the  Hewlett-Packard  9000  computer. 
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10  !  COMPUTE  COVflRIflKCE  MATRICES  AMD  EIGENVALUES  FOR  WHITE 
20  !  NOISE  INPUT;  STORE  EIGENVALUES  IN  "EIG”  IN  LINE  1720 
30  Ka=-10  !  INPUT  SIGNAL  START 

40  Kb=5  !  INPUT  SIGNAL  END;  Kb  >=  Ka 

50  Kc=0  !  ACCUMULATOR  START 

60  Kd=5  !  ACCUMULATOR  END;  Kd  >=  Kc,Ka 

?0  Delta=.2  !  TIME  SAMPLING  INCREMENT  (SECONDS 

80  J=3  I  NUMBER  OF  SIGNAL  COMPONENTS 

90  DATA  11., -48., 75,  t  SIGNAL  SCALINGS  <WE  SET  SUM  =  1> 

100  DATA  1  .,.  125,  .  066667  !  SIGNAL  TIME  CONSTANTS  (SECONDS.) 

110  N=4  !  HUMBER  OF  FILTER  COMPONENTS 

120  DATA  l.,2.,3.,4.  !  FILTER  SCALINGS  (ARBITRARY) 

130  DATA  .3,. 5,. 7,. 9  I  FILTER  TIME  CONSTANTS  (SECONDS) 

140  Nr=20  I  NUMBER  OF  SIGNAL-TO-NOISE  RATIOS 

150  DATA  0,1,2,3,4,5,6,7,8,9,10  !  INPUT  SNRS  IN  DB 

160  DATA  11,12,13,14,15,16,17,18,19 

170  IF  (Ka<=Kb)  AND  (Kc<=Kd)  AND  (Ka<=Kd)  THEN  200 

180  PRINT  "PROBLEM  WITH  PARAMETERS" 

190  STOP 

200  DIM  A1 pha( 10) , A( 10) , Psi < 10) , C( 10) , R(20) , N( 100) 

210  DIM  Pc(10),Ec(10),Pc<.)(10),Cn(200),Ea(10),E(0,  1000) 

220  DIM  Gca(100),Gcda(100),Gcb(100),Gbda( 100),Cs(5000) 

230  DIM  Cw(0, 10000), Ca<0, 10000), D(100),Eig(0, 2000) 

240  DOUBLE  Ka, Kb , Kc , Kd , J , N , Nr , Ns , Ks , Kdc , Ko  !  INTEGERS 
250  DOUBLE  L 1 , L2 , L , L 1 1 ,  Js , Qs  ,  K i  ,  Ks 1 ,  K 1 ,  Ks 2  ,  K2  ,  I 

260  REDIM  Alphad!  J),A(1  j  J) 

270  READ  Alpha(*) , A<*) 

280  MAT  A1pha=Alphar(SUM( Alpha)) 

290  MAT  A=(DeUa)/A 

300  REDIM  Psi ( 1 sN) , C< 1 1 N) 

310  READ  Psi (»),C(*) 

320  MAT  C  =  (DeUa)-'C 

330  REDIM  R(l:Mr) 

340  READ  R(*) 

350  RFDIM  M<Kc;Kd) 

360  CALL  We i ght s ( Kc , Kd , W < * ) ) 

370  MAT  W=SQR(W) 

380  REDIM  Pc ( 1 : N) , Ec ( 1 ; N)  !  FILTER  OUTPUT  NOISE  COVARIANCE 
390  FOR  Ns=l  TO  N 

400  C=C(N5) 

410  Pc (Ns)=Psi (Ns)*C 

420  Ec (N5)=EXP(-C) 

430  NEXT  Ns 

440  RED  Ml  Pcm(1:N) 

450  FOR  Ns=l  TO  N 

460  E=Ec(Ns) 

470  S=0. 

480  FOR  Ks=l  TO  N 

490  S-C  +  rc <Ks)/( 1 . -E*Ec (Ks)  ) 

500  NEXT  Ks 

510  Pcc<(Ns)=Pc  (Ns)*S 

520  NEXT  Ns 

530  Kdc=Kd-Kc 

540  REDIM  Cn(0;Kdc) 

550  HAT  Cn=(0.) 

560  FOR  Ns=l  TO  N 

570  Pcv)  =  Pco(Ns) 

580  E=Ec(Ns) 
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590  P=l. 

600  Cn  <  0 )  =Cn  <  0 ) +Pc 

610  FOR  Ks=l  TO  Kdc 

620  P=P*E 

630  Cn<Ks)=Cn<K£)+Pcv*P 

640  NEXT  Ks 

650  HEXT  Ns 

660  FOR  Js=l  TO  J  !  FILTER  OUTPUT  SIGNAL  COVARIANCE 

670  A< Js)=EXP<R< Js) ) 

680  NEXT  Js 

690  Ko=MAX<Ka, Kc ) 

700  REDIM  Ea<l;N),E<l:H,KojKd) 

710  FOR  Ns=l  TO  N 

720  E=Ec<Ns) 

730  Ea<Ns)=l./E 

740  E<Ns,Ko)=EXP<-C<Ns)*Ko) 

750  FOR  Ks=Ko+l  TO  Kd 

760  E<Ns,Ks)=E<Ns,Ks-l )*E 

770  NEXT  Ks 

780  NEXT  Ns 

790  Ll=MIN<Ko,Kb)+l 

800  L2=MIN<Kd,Kb)+l 

810  REDIM  Gca(Ll:L2),Gcda<:Ll:L2),Gcb<:Ll:L2>,Gbda<Ll:L2> 

820  L=Kd-Ko+l 

830  REDIM  Cs< 1 ; L*<L+1 )/2) 

840  Lll=Ll+l 

850  FOR  Js=l  TO  J 

860  Al=A1pha<Js) 

870  A=A<Js) 

880  A1*A*A-1. 

890  FOR  Ns=l  TO  N 

900  C=Ea<Ns) 

910  Ca=C»A 

920  Cda=C/A 

930  Cal=l.-Ca 

940  Ac=A-C 

950  Ae=Al*Pc<Ns) 

960  Gc  a=Ca''Ka/Cal 

970  GcaCLl  )=Ca-Ll/'Cal 

980  Gcda<Ll )=Cda^Ll/< 1 . -Cda) 

990  FOR  Ks=Lll  TO  L2 

1000  Gca<Ks)=Ca*Gca<Ks-l ) 

1010  Gcda<Ks)=Cda*Gcda<Ks-l ) 

1020  NEXT  Ks 

1030  FOR  Qs=l  TO  N 

1040  B=Ea<Qs> 

1050  Aee  =  Ae*Pc < Qs ) 

1060  Ba=>B«A 

1070  Bda=B/A 

1080  Cb=C*E 

1090  Bal=l.-Ba 

1100  Cbl=l.-Cb 

1110  Gba=Ba''Ka/Bal 

1120  Rb=A-B 

1130  F=B*Al/<Ab*Bal > 

1  140  Sl=Cb''Ka/CbHKAl+Cbl  )/<Ab*Ac  ) 

1150  GcbCLl )=Cb^Ll/Cbl 

1160  Gbda<Ll  )=Bda^Ll/"<  1 . -Bda) 
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1170  FOR  Ks=Lll  TO  L2 

1180  Gcb<Ks)=Cb*Gcb<Kt-l > 

li;?  Gbdac:Ks)=Bda*Gbda<Ks-l  ) 

1200  NEXT  Ks 

12  0  K’=0 

1220  FOR  K5l=Ko  TO  Kd 

1230  Kl=MIN<Ksl , Kb)+1 

1240  S2  =  Sl+Gcb(Kl  )*F-Gba*G.:da<Kl  ) 

1250  S3  =  Gca<:Kl  )-Gca 

1260  E  =  Ree*E<*!s,  Ksl ) 

127 ,  FOR  Ks2=Ksl  TO  Kd 

1280  Ki=Ki+l 

1290  K2  =  MIN<i<s2,  Kb)  +  1 

1300  S=S2+Gbda<K2)*S3  ‘  SUM  S 

1310  CsCKi )=Cs<Ki ) +E*E < Qs , Ks2 ) *3 

K 20  NEXT  Ks2 

1330  NEXT  Ksl 

1340  NEXT  Qs 

1350  NEXT  Ns 

1360  NEXT  Js 

i3P0  L  =  Kd-Kc  +  l  I  TOTAL  NEIGHTED  COVFIRIftNCE  MATRIX 

1380  REDIM  Cw< 1 : L, 1 ; L) , Ca< 1 ! L, 1 ; L) , D< 1 : L) , El g<0; Nr , 1 : L) 

1390  Lll=^Kc-l 

1 400  R=0 • 

1410  FOR  1=0  TO  Nr 

1420  K1=0 

1430  IF  1=0  THEN  1450 

1440  R=10. ^<R< D*. 1  )  !  INPUT  POWER  S I GNAL-TO-NO I SE  RATIO 

1450  FOR  Ksl=Kc  TO  Kd 

1460  W=W<Ksl)  (  WEIGHTS  (u(k>> 

1470  FOR  Ks2=Ksl  TO  Kd 

1480  IF  Ksl>=Ko  THEN  1510 

1490  Cs=0. 

1500  GOTO  1530 

15*0  Ki=Ki+l 

1520  Cs=Cs<Ki) 

1530  Pr=W*W<Ks2)*<Cn<Ks2-Ksl )+R*C£) 

1540  Ll=Ksl-Lll 

1550  L2=K32-L11 

1560  Cw<Ll , L2i=Cw(L2, LI )=Pr 

1570  NEXT  Ks2 

1580  NEXT  Ksl 

1590  CALL  Swd<L,L,Cw<*),Ca<»>,D<*))  !  EIGENVALUES 

1600  MAT  SORT  DC*)  DES 

1610  IF  D<L)>0.  THEN  1640 

1620  PRINT  "PROBLEM:  SOME  NON-POSITIVE  EIGENVALUES" 

1630  PAUSE 

1640  PRINT  "I  ="jl5"  CONDITION  NUMBER  ="jDCl)/DCL> 

1650  PRINT  DC*) 

1660  PRINT 

1670  FOR  K£=1  to  L 

1680  Ei gC  I  , Ks)=DCKs)  !  STORE  EIGENVALUES 

1690  NEXT  Ks 

1700  NEXT  I 

1710  MASS  STORAGE  IS  ":CS80,7" 

1720  ASSIGN  #l  TO  "EIG" 

1730  PRINT  #1 ; Nr , Kc , Kd, Ei gC*) 

1740  ASSIGN  #1  TO  * 

1750  END 

1760  ! 


62 


TR  8913 


1770  SUB  Weight sCDOUBLE  Kc,Kd,REflL  g<*>> 

1780  DOUBLE  Ks  I  INTEGER 

1790  F=l.  !  DECAY  FACTOR  (DIMENSIONLESS) 

1800  W<Kd)=l. 

1810  FOR  Ks=Kd-l  TO  Kc  STEP  -1 

1820  W(Ks)=W<Ks+l)*F  !  EXPONENTIAL  WEIGHTS 

1830  NEXT  Ks 

1840  SUBEND 

1850  ! 

1880  SUB  SMdCDOUBLE  M,N,REAL  A< * ) , V < * ) , W ( * ) ) 

1870  ALLOCATE  R^^KliN)  !  NUMERICAL  RECIPES,  PAGES  60-64 

1880  IF  M>=N  THEN  1910  !  A<*)  IS  OVER-WRITTEN 

1890  PRINT  '‘M<N  IS  DISALLOWED" 

1900  PAUSE 

1910  DOUBLE  I, J,K,L, Its,Nm, Jj  I  INTEGERS  (NOT  DOUBLE  PRECISION 
1920  G  =  Sc  al  e  =  Anortri  =  0 . 

1930  FOR  1=1  TO  N 

1940  L=I+1 

1950  Ryl(I)=Scale*G 

I960  G=S=Sca1e=0. 

1970  IF  I>M  THEN  2250 

1980  FOR  K=I  TO  M 

1990  Sca1e  =  Sca1e  +  ABS(P(K,  D) 

2000  NEXT  K 

2010  IF  Scale=0.  THEN  2250 

2020  FOR  K=I  TO  M 

2030  Aa=A(K, I)=A(K, I)/Scale 

2040  S=S+Aa*Aa 

2050  NEXT  K 

2060  F=A(I,I) 

2070  G=-SQR(S) 

2080  IF  F<0.  THEN  G=-G 

2090  H=F*G-S 

2100  A(I,I)=F-G 

2110  IF  I=N  THEN  2220 

2120  FOR  J=L  TO  N 

2130  S=0. 

2140  FOR  K=I  TO  M 

2150  S=S+A(K, I >*A(K, J) 

2160  NEXT  K 

2170  F=S/H 

2180  FOR  K=I  TO  M 

2190  A(K, J)=A(K, J)+F*A(K, I  ) 

2200  NEXT  K 

2210  NEXT  J 

2220  FOR  K=I  TO  M 

2230  A(K, I )=A(K, I)*Scale 

2240  NEXT  K 

2250  W(I)=Scale*G 

2260  G=S=Scale=0. 

2270  IF  (I>M)  OR  (I=N)  THEN  2570 

2280  FOR  K=L  TO  N 

2290  Seal e=Sca1 e+RBS(A( I , K) ) 

2300  NEXT  K 

2310  IF  Scale=0.  THEN  2570 
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2320  FOR  K=L  TO  N 

2330  fia=fl<I,K>=fl<I,K)/Scale 

2340  S=S+fla*fla 

2350  NEXT  K 

2360  F=R<I,L> 

2370  G=-SQR<S) 

2380  IF  F<0.  THEN  G=-G 

2390  H=F*G-S 

2400  fl(I,L)=F-G 

2410  FOR  K=L  TO  N 

2420  RmI  <:K)=FI<  I , 

2430  NEXT  K 

2440  IF  I=M  THEN  2540 

2450  FOR  J=L  TO  M 

2460  S=0. 

2470  FOR  K=L  TO  N 

2480  S=S+fl< J, K)*fl< I , K) 

2490  NEXT  K 

2500  FOR  K=L  TO  N 

2510  J,  K)=fl<  J,  K)+S*R>j1  <K) 

2520  NEXT  K 

2530  NEXT  J 

2540  FOR  K=L  TO  H 

2550  fl( I , K)=fl< I , K)*Scal e 

2560  NEXT  K 

2570  flnor  rn  =  MFIXCfinorfn,  flBG^lK  I  )  )+ftBS<R<^l  <  I )  )  > 

2580  NEXT  I 

2590  FOR  I=N  TO  1  STEP  -1 

2600  IF  I>=N  THEN  2770 

2610  IF  G=0.  THEN  2740 

2620  FOR  J=L  TO  N 

2630  V<J, I)=R<I, J)/fl<I,L)/G 

2640  NEXT  J 

2650  FOR  J=L  TO  N 

2660  S=0. 

2670  FOR  K=L  TO  N 

2680  S=S+fl< I , K)*V<K, J) 

2690  NEXT  K 

2700  FOR  K=L  TO  N 

2710  V<K, J)=V<K, J)+S»V<K, I  ) 

2720  NEXT  K 

2730  NEXT  J 

2740  FOR  J=L  TO  N 

2750  V<  I  ,  .1)=^;  J,  I  )=0. 

2760  NEXT  J 

2770  V<:i,I)  =  l. 

2780  G  =  Rm1<I) 

2790  L=I 

2800  NEXT  I 

2810  FOR  I=N  TO  1  STEP  -1 

2820  L=I+1 

2830  G=W(I) 

2840  IF  I>=N  THEN  2880 

2850  FOR  J=L  TO  N 

2860  fi<I,.J)=0. 

2870  NEXT  J 
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2880  IF  G=0.  THEN  3050 

2890  G=l./G 

2900  IF  I=N  THEN  3010 

2910  FOR  J=L  TO  N 

2920  8=0. 

2930  FOR  K=L  TO  M 

2940  S  =  S  +  fl<K,  I  J) 

2950  NEXT  K 

2960  F  =  S/'fl<  I ,  I  )*G 

2970  FOR  K=I  TO  M 

2980  H<K, J)=fl<K, J>+F*fl<K, I  ) 

2990  NEXT  K 

3000  NEXT  J 

3010  FOR  J=I  TO  M 

3020  R<  J,  I  >=fl<  J,  n*G 

3030  NEXT  J 

3040  GOTO  3080 

3050  FOR  J=I  TO  M 

3060  fl<.J,I)=0. 

3070  NEXT  J 

3080  fl< I , I )=fi< I , I)+l . 

3090  NEXT  I 

3100  FOR  K=N  TO  1  STEP  -1 
3110  FOR  Its=l  TO  30 

3120  FOR  L=K  TO  1  STEP  -1 
3130  Nni  =  L-l 

3140  IF  <flBS<Rvl  <L)  )+flnorffi)=flnortn  THEN  3360 
3150  IF  <flBS<W<Nm)  )  +  flnorfn)  =  flnorm  THEN  3170 
3160  NEXT  L 

3170  C=0. 

3180  S=l. 

3190  FOR  I=L  TO  K 

3200  F  =  S*Rvl<n 

3210  Rvl < I)=C*Rvl < I ) 

3220  IF  < flBS < F  ) +finorrrt  )  =flriorm  THEN  3360 
3230  G=W<I) 

3240  H=SQR<F*F+G*G) 

3250  W<I)=H 

3260  H=l./H 

3270  C=G«H 

3280  S=-F*H 

3290  FOR  J=1  TO  M 

3300  Y=fi<J,Nm) 

3310  2=fl<J,n 

3320  fl< J, Nm)=Y*C+Z*S 

3330  FIC.J,  I  )=-Y*S  +  2*C 

3340  NEXT  J 

3350  NEXT  I 

3360  2=W<K) 

3370  IF  LOK  THEN  3440 

3380  IF  2>=0.  THEN  3430 

3390  W':K>=-2 

3400  FOR  J=1  TO  N 

3410  V(  J,K>=-V(:  J,K) 

3420  NEXT  J 

3430  GOTO  3970 
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3440  IF  Its<30  THEN  3470 

3450  PRINT  "H0  CONVERGENCE  IN  30  ITERATIONS" 

3460  PAUSE 

3470  X=W<L) 

3480  Nro=K-l 

3490  Y=W<Nm) 

3500  G  =  R^.M<Nm) 

3510  H=Rul<K) 

3520  F=<  <Y-2)*<Y  +  Z)  +  <:G-H)*<G  +  H))/<2.  ♦H*Y) 

3530  G  =  SQR<F*F+l  .  ) 

3540  Aa=ABS<G) 

3550  IF  F<0.  THEN  Aa=-Aa 

3560  F=';  <X-Z)*CX  +  2)+H*<  <Y/<F  +  Aa)  )-H)  )/X 

3570  C=S=1. 

3580  FOR  J=L  TO  Nm 

3590  I=J+1 

3600  G=Rm1<I) 

3610  Y=W<I) 

3620  H=S*G 

3630  G=C*G 

3640  Z=SQR<F*F+H*H) 

3650  R^^1<J)=2 

3660  C  =  F-^Z 

3670  S=H/Z 

3680  F=X*C+G*S 

3690  G=-X*S+G»C 

3700  H=Y*S 

3710  Y=Y*C 

3720  FOR  .J,j  =  l  TO  N 

3730  X=V<Jj,J> 

3740  2=V<Jj,I) 

3750  V( J j , J)=X*C+Z*S 

3760  V<.J  j  ,  I)=-X»S  +  Z*C 

3770  NEXT  Jj 

3780  Z=SQR<F»F+H*H) 

3790  W(J)=Z 

3800  IF  Z=0.  THEN  3840 

3810  Z=l./Z 

3820  C=F*Z 

3830  S=H*Z 

3840  F=C*G+S«Y 

3850  X=-S»G+C*Y 

3860  FOR  Jj=l  TO  M 

3870  Y=A<Jj,J) 

3880  Z  =  A(J.j,I> 

3890  A< J j , J)=Y*C+Z*S 

3900  AC J j , I )=-Y*S+Z*C 

3910  NEXT  Jj 

3920  NEXT  J 

3930  Rm1<L)=0. 

3940  Rvl<K)=F 

3950  W<K)=X 

3960  NEXT  Its 

3970  NEXT  K 

3980  SUBEND 
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10 

!  PLOT  ROC  FROM  EIGENVALUES 

"EIG",  LINES  120  8.  1060 

20 

De 1 X i = .  08  1 

INITIAL  INCREMENT  ON  CHAR.  FUNCTION 

30 

M=128  ! 

INITIAL  SIZE  OF  FFT 

40 

Bs=0.  I 

SHIFT  b 

50 

Tol0=l.E-36  ! 

TOLERANCE  FOR  FALSE  ALARM  PROBABILIT 

60 

Toll=l.E-32  ! 

TOLERANCE  FOR  DETECTION  PROBABILITIE 

70 

DOUBLE  M, Nr, Kc , Kd, Kdc , Ms, 

I,Ks,Ns,  Im,Mivi  !  INTEGERS 

80 

DIM  E(0: 20, 1 ;  100) , D< 1  I  100) , Co£<0; 4096)  !  20  NONZERO  SNRs 

90 

DIM  X<0i  16383), Y<0: 16383) 

, Pr<0: 20, 0: 127) 

100 

PRINT  "Delxi  = " ; De 1 x i ,  "M 

=“;M,"b  ="jBs 

110 

MASS  STORAGE  IS  "rCSSO,?" 

120 

ASSIGN  #1  TO  "EIG" 

130 

READ  #l;Nr,Kc,Kd  ! 

NO.  OF  SNRs,  START  8:  END  OF  SUMMER 

140 

Kdc=Kd-Kc+l 

150 

REDIM  E<0;Nr, l!Kdc),D<l:Kdc> 

160 

READ  #l5E<*)  1 

EIGENVALUES 

170 

GINIT 

180 

PLOTTER  IS  "GRAPHICS" 

190 

GRAPHICS  ON 

200 

REDIM  Co£<0; Mx4) , X<0; M-1  ) 

, Y<0: M-1 ) 

210 

A=2. »PIxM 

220 

FOR  Ms=0  TO  M/4 

230 

Cos<Ms)=COS<A*Ms)  ! 

INITIAL  QUARTER-COSINE  TABLE 

240 

NEXT  Ms 

250 

Tol=Tol0 

260 

FOR  1=0  TO  Nr  ! 

Nr  NONZERO  S I GNAL-TO-NO I SE  RATIOS 

270 

IF  I >0  THEN  Tol =To1 1 

280 

FOR  Ks=l  TO  Kdc 

290 

D<Ks)=E< I , Ks)  ! 

COPY  EIGENVALUES 

300 

NEXT  Ks 

310 

Mux=SUM<D)  ! 

MEAN  OF  RANDOM  VARIABLE  x 

320 

R  =  0.  ! 

ARGUMENT  OF  SQUARE  ROOT 

330 

P=1 .  1 

POLARITY  INDICATOR 

340 

Muy=Mux+Bs  ! 

MEAN  OF  y  =  X  +  b 

350 

MAT  X=(0. ) 

360 

MAT  Y=<0. ) 

370 

Y<0)=. 5*De1 xi *Nuv 

380 

Ns  =  0 

390 

Ns=Ns+l  1 

LOOP  ON  xi 

400 

Xi=Delxi*Ns  1 

ARGUMENT  xi  OF  CHAR.  FN. 

410 

Txi=Xi*2.  1 

CALCULATION  OF 

420 

Pr=l.  1 

CHARACTERISTIC 

430 

Pi=0.  ! 

FUNCTION  fy<xi) 

440 

FOR  Ks=l  TO  Kdc 

450 

Te  =  Tx i *0  <  Ks ) 

460 

Pe=Pr+Te»Pi 

470 

P i =P i -Te*Pr 

480 

Pr  =  Pe 

490 

NEXT  Ks 

500 

CALL  Sqr < Pr , P i , Sr , S i  ) 

510 

De=Sr»Sr+Si ♦Si 

520 

F  yr  =  Sr /De 

530 

F  y i =-S i /De 

540 

Ro  =  R 

550 

R=ATN(Fy i /Fyr ) 

560 

IF  RBS<R-Ro) >1 . 6  THEN  P=- 

■P 

570 

IF  P>0.  THEN  600 
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580  Fyr =-Fyr 

590  Fyi=-Fyi 

600  Ms=Ns  MODULO  M  !  COLLflPSIHG 

610  fl=Fyr/Ns 

620  B  =  Fy1/'Ns 

630  X<Ms)=X<Mi)+fl 

640  Y<Ms)=Y<Ms)+B 

650  IF  fl*fl+B*B<Tol  THEN  670 

660  GOTO  390 

670  CALL  Fftl4<M,Cos<*>,X<*),Y<»)> 

680  GCLEAR 

690  WINDOW  0,M,-18,0 

700  LINE  TYPE  3 

710  GRID  M/8,  1 

720  LINE  TYPE  1 

730  FOR  Ks=0  TO  M-1 

740  T=Y<Ks)/PI-Ks/M 

750  Y<Ks)=Pr=.  5  +  T  !  EXCEEDANCE  PROBABILITY  IN  Y<*;:' 

760  IF  Pr>=l.E-16  THEN  Y  =  LGT  <  Pr  )  !  Y  <  Ks  )  =PROB  <  x  >2  PI  Ki-'CM  De  1  x  i  > -Bs  > 

770  IF  Pr<=-1.E-16  THEN  Y=-32. -LGT<-Pr ) 

780  IF  ABS<Pr')<  l.E-16  THEN  Y  =  -16. 

790  PLOT  Ks,Y 

800  NEXT  Ks 

810  PENUP 

820  IF  I=Mr  THEN  940 

830  PRINT  Ij 

840  INPUT  "SCALE  FFT  SIZE  BY  1 , 2 , 4 , 8 ,  .  .  .  :  "  ,  I rn 

850  Mm=MIN<M*Im, 16384) 

860  IF  Mm=N  THEN  940 

870  De 1 X i =De 1 X i *M/Nm 

880  M  =  Mm 

890  REDIM  Cos<0:M/4),X<0:M-l),Y<0:M-l) 

900  A=2.*PI/M 

910  FOR  M£=0  to  M/4 

920  Cos<Ms)=COS<R*Ms)  !  QUARTER-COSINE  TABLE;  <=  4096 

930  NEXT  Ms 

940  FOR  Ks=0  TO  127 

950  Pr<I,K3)=Y<Ks)  !  STORE  FOR  ROC  PLOT 

960  NEXT  Ks 

970  NEXT  I 

980  FOR  Ks=0  TO  127 

990  IF  Pr<0,Ks)<.l  THEN  1010 

1000  NEXT  Ks 

1010  Ns  =  Ks-l 

1020  FOR  Ks=Ms  TO  127 

1030  IF  Pr<0,Ks)< lE-10  THEN  1050 

1040  NEXT  Ks 

1050  Ns=Ks 

1060  ASSIGN  #1  TO  "ROC" 

1070  PRINT  #1; Ms, Ns, Nr 

1080  PRINT  #l5Pr<*) 

1090  ASSIGN  #1  TO  * 

1100  DIM  Pf < 1 : 14) , Pd< 1 j 18) 

1110  DATA  lE-10, lE-9, lE-8, lE-7, lE-6, lE-5, lE-4 
1120  DATA  . 001 , . 002, . 005, . 01 , . 02, .05, . I 

1130  READ  Pf<*) 
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1140  DflTfi  .01, .02, .05, . 1, .2, .3, .4, .5, .6, .7 
1150  DATA  .8, .9, .95, .98, .99, .995, .998, .999 
1160  READ  Pd<*) 

1170  FOR  1=1  TO  14 

1180  Pf<I)=FNInMphi 

1190  NEXT  I 

1200  FOR  1=1  TO  18 

1210  Pd<I>=FNInvphi <Pd<I)) 

1220  NEXT  I 

1230  GCLEAR 

1240  Xl=Pf<l) 

1250  X2=Pf<14) 

1260  Yl=Pd(l) 

1270  Y2=Pd<18) 

1280  WINDOW  X1,X2,Y1,Y2 

1290  LINE  TYPE  3 

1300  FOR  1=1  TO  14 

1310  MOVE  Pf<I),Yl 

1320  DRAW  Pf<n,Y2 

1330  NEXT  I 

1340  FOR  1=1  TO  18 

1350  MOVE  Xl,Pd<I) 

1360  DRAW  X2,Pd<I) 

1370  NEXT  I 

1380  PENUP 

1390  LINE  TYPE  1 

1400  FOR  Ks=Ms  TO  Ns 

1410  Pr<0,Ks)=FNInvphi <Pr<0,Ks))  !  FALSE  ALARM  PROBABILITY 

1420  NEXT  Ks 

1430  FOR  1=1  TO  Nr 

1440  FOR  Ks«Ms  TO  Ns 

1450  X=Pr<0,Ks) 

1460  Pr=Pr<I,Ks) 

1470  IF  Pr>.9999  THEN  1510 

1 480  Pr  < I , Ks ) =Y  =  FN Inwph i <Pr ) 

1490  PLOT  X,Y 

1500  IF  Pr<.01  THEN  1520 

1510  NEXT  Ks 

1520  PENUP 

1530  NEXT  I 

1540  PAUSE 

1550  END 

1560  I 

1570  DEF  FNArg<X,Y)  !  PRINCIPAL  ARGCZ) 

1580  IF  X=0.  THEN  RETURN  .5*PI*SGN<Y) 

1590  A=ATN<Y/X) 

1600  IF  X>0.  THEN  RETURN  A 

1610  IF  Y<0.  THEN  RETURN  A-PI 

1620  RETURN  A+PI 

1630  FNEND 

1640  ) 

1650  SUB  Sqr<X,Y,U,V)  I  PRINCIPAL  SQR<2) 

1660  F=SQR<SQR<X*X+Y*Y) ) 

1670  T=. 5*FNArg<X, Y) 

1680  U=F*COS<T) 

1690  V=F*SIN<T) 

1700  SUBEND 

1710  ! 
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1720  DEF  FHInMphi<X)  !  AMS  55,  26.2.23 

1730  IF  X=.5  THEN  RETURN  0. 

1740  P=MIN(X, l.-X) 

1750  P  =  riflX':P,  l.E-20> 

1760  T=-LOG<P) 

1770  T=SQR<T+T) 

1780  P=l.+T*<1.432788+T»<. 189269+T*. 001 308> ) 

1790  P=T-<2. 515517+T*< . 802853+T*. 010328) )/P 

1800  IF  X<.5  THEN  P=-P 

1810  RETURN  P 

1820  FNEND 

1830  ! 

1840  SUB  Fftl4<D0UBLE  N,REflL  Cos  <  *  )  ,  X  <  *  )  ,  Y  <  *  )  )  !  .J<  =2  ■  1  4=  1 6384  ;  0  SUBS 
1850  DOUBLE  Log2n ,  N 1  ,  N2 ,  N3 ,  N4 ,  J ,  K  !  INTEGERS  <  2''31  =  2,147,483,648 

I860  DOUBLE  I  1 ,  I  2 ,  I  3 , I  4 ,  I  5 , 1 6 , 1 7 , 1 8 , 1 9 , I  1 0 , I  1 1 ,  I  1 2 , I  1 3 ,  I  1 4 , L ^ 0 :  1 3 ) 

1870  IF  N=1  THEN  SUBEXIT 

1880  IF  N>2  THEN  1960 

1890  fi  =  X<0)+XU) 

1900  X< 1 )=X<0)-X( 1  ) 

1910  X<0)=fl 

1920  fl  =  Y<:0)+Y<l) 

1930  Y< 1 )=) (0)-Y<  1  ) 

1940  Y<0)=fl 

1950  SUBEXIT 

1960  fi  =  L0G<N)/L0G<2.  ) 

1970  Log2n=fl 

1980  IF  flBS<fl-Log2n)< 1 . E-8  THEN  2010 

1990  PRINT  "N  ="jN}"IS  NOT  A  POWER  OF  2;  DISALLOWED." 

2000  PAUSE 

2010  Hl=N/4 

2020  N2=N1+1 

2030  N3=H2+1 

2040  N4=N3+N1 

2050  FOR  11=1  TO  Log2n 

2060  I2  =  2-':Log2n-1 1 ) 

2070  13=2*12 

2080  I4=N/I3 

2090  FOR  15=1  TO  12 

2100  I6=< 15-1  )*I4+1 

2110  IF  I6<=N2  THEN  2150 

2120  Rl=-Cos<N4-I6-l ) 

2130  A2=-Cos< I6-N1-1 ) 

2140  GOTO  2170 

2150  Al=Cos':  16-1  ) 

2160  A2=-Cos<N3-I6-l ) 

2170  FOR  17=0  TO  N-I3  STEP  13 

2180  18=17+15-1 

2190  19=18+12 

2200  T1=X<I8) 

2210  T2=X<I9) 

2220  T3=Y<I8) 

2230  T4=y<19) 
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2240 

2250 

2260 

2270 

2280 

2290 

2300 

2310 

2320 

2330 

2340 

2350 

2360 

2370 

2380 

2390 

2400 

2410 

2420 

2430 

2440 

2450 

2460 

2470 

2480 

2490 

2500 

2510 

2520 

2530 

2540 

2550 

2560 

2570 

2580 

2590 

2600 

2610 

2620 

2630 

2640 

2650 

2660 

2670 

2680 

2690 

2700 

2710 

2720 

2730 

2740 

2750 

2760 

2770 


R3=T1-T2 

84=13-74 

X< I8)=T1+T2 

Y< I8)=T3+T4 

X<  I9>=Rl*fl3-fl2*Fl4 

Y<  I9)=fil»fi4  +  R2*fl3 

NEXT  17 

NEXT  15 

NEXT  II 

1 1 =Log2n+ 1 

FOR  12=1  TO  14 

L<  12-0  =  1 

IF  I2>Log2n  THEN  2380 
L(I2-0=2''<I1-I2) 

NEXT  12 
K  =  0 

FOR  11=1  TO  L<13) 

FOR  12=11  TO  L<12)  STEP  L<13) 
FOR  13=12  TO  L<11)  STEP  L<12) 
FOR  14=13  TO  L<i0)  STEP  L<11) 
FOR  15=14  TO  L<9>  STEP  L<10> 
FOR  16=15  TO  L<8)  STEP  L<9) 
FOR  17=16  TO  L<7)  STEP  L<8) 
FOR  18=17  TO  L<6)  STEP  L<7> 
FOR  19=18  TO  L<5)  STEP  L<6) 
FOR  110=19  TO  L<4)  STEP  L<5) 
FOR  111=110  TO  L<3)  STEP  L<4) 
FOR  112=111  TO  L<2)  STEP  L<3) 
FOR  113=112  TO  LCD  STEP  L<2) 
FOR  114=113  TO  L<0)  STEP  L<1) 
J=I14-1 

IF  K>J  THEN  2620 
R=X<K> 

X<K)=X< J) 

X< J)=R 
R=Y<K> 

Y<K)=Y< J) 

Y< J)=R 
K  =  K+1 
NEXT  114 
NEXT  113 
NEXT  112 
NEXT  Ill 
NEXT  110 
NEXT  19 
NEXT  18 
NEXT  17 
NEXT  16 
NEXT  15 
NEXT  14 
NEXT  13 
NEXT  12 
NEXT  II 
SUBEND 
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10 

!  SIMULATION  FOR  ONE-POLE 

PROCESSES 

20 

A5=.85  ! 

INPUT  SIGNAL 

PARAMETER 

30 

Bs=.77  ! 

INPUT  NOISE 

PARAMETER: 

40 

Cs=1.13  ! 

FILTER  PARAMETER 

50 

R=.  7  ! 

INPUT  POWER  : 

SNR 

60 

Ka=4 

70 

Kb=l  1 

80 

Kc=2 

90 

Kd=16 

100 

Bi ns=1000 

110 

Hum=lE6  I 

NUMBER  OF  TRIALS 

120 

Deltaz=.25  ! 

BIN  SIZE  FOR 

OUTPUT  Z 

130 

DOUBLE  Ka,Kb,Kc,Kd,Bins 

, Num, It , K, M 

140 

DIM  E<  100) , 3< 100) , Yn<100), YsC 100>, T( 

1  000  ) 

150 

Ea=EXP(-As) 

160 

Eb=EXP<-Bs)  ! 

WHITE  NOISE: 

FOR  Bs-uif 

170 

Ec=EXP<-Cs> 

180 

Fa=SQR(R*< 1 . -Ea»Ea) ) 

190 

Fb=SQR< 1 . -Eb*Eb) 

200 

REDIM  E<0jKd-Ka) 

210 

E<0)=Cs 

220 

FOR  K=1  TO  Kd-Ka 

230 

E<K)=E<K-1 )*Ec 

240 

NEXT  K 

250 

REDIM  S<-30:Kd),Yn<Kc-l 

; Kd), YsCKa: Kd ) , T < 1 : B i ns  ) 

260 

FOR  It=l  TO  Num 

270 

S<-30>=Nk*Yk=0. 

280 

FOR  K=-29  TO  Kc-1  ! 

ALLOW  STEADY 

STATE 

290 

R1=2.»RND-1.  ! 

BOX-MULLER 

300 

R2  =  2. *RND-1 . 

310 

R3=R1*R1+R2*R2 

320 

IF  R3>1.  THEN  290 

330 

R3=SQR<-2. *L0G<R3)/R3) 

340 

R1=R1*R3 

350 

R2=R2*R3 

360 

S<K)=Fa*Rl+Ea*S<K-l )  ! 

FILTER  INPUT 

SIGNAL 

370 

Nk=Fb*R2+Eb*Nk  ! 

FILTER  INPUT 

NOISE 

380 

Yk=Cs*Hk+Ec*Yk  ! 

FILTER  OUTPUT  NOISE 

390 

NEXT  K 

400 

Yn<Kc-l )=Yk 

410 

FOR  K=Kc  TO  Kd 

420 

Rl=2. *RND-1 . 

430 

R2=2. *RND-1 . 

440 

R3  =  Ri*Rl-t-R2*R2 

450 

IF  R3>1.  THEN  420 

460 

R3=SQR<-2. *L0G<R3)/R3) 

470 

R1=R1*R3 

480 

R2=R2*R3 

490 

M  =  K-1 
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500 

S<K>=Fii*Rl+Ea*S<M)  1 

FILTER 

INPUT 

510 

Nk=Fb*R2+Eb*Nk  ! 

FILTER 

INPUT 

520 

Yn<K)=Cs*Nk+Ec*Yn<M)  1 

FILTER 

OUTPUT 

530 

NEXT  K 

540 

FOR  K=Ka  TO  Kd 

550 

S  =  0. 

560 

FOR  M=Ka  TO  MIN<K,Kb) 

570 

S  =  S  +  E<K-M)*S':M)  I 

FILTER 

OUTPUT 

580 

NEXT  M 

590 

Ys<K)=S 

600 

NEXT  K 

610 

Z  =  0. 

620 

FOR  K=Kc  TO  Ka-1 

630 

T=Yn<K) 

640 

Z  =  Z  +  T*T 

650 

NEXT  K 

660 

FOR  K=Ka  TO  Kd 

670  T=Ys<K)+Yn<K) 

680  Z=2+T*T 

696  HEXT  K 

700  K=IHT<Z/Del taz)+l 

710  K=MIN<K,Bins) 

720  T(K)=T<K)+1. 

730  NEXT  It 

740  FOR  K=Bin£-l  TO  1  STEP  -1 
750  T<K)=T<K)  +  T<K+1  ) 

760  NEXT  K 

770  MAT  TsT/dluifi) 

780  GINIT 

790  PLOTTER  IS  "GRAPHICS" 

800  GRAPHICS  ON 

810  WINDOW  0, Bins, -5,0 

820  GRID  B i  ns/'4 ,  1 

830  FOR  K=1  TO  Bins 

840  T=T<K) 

850  IF  T>0.  THEN  870 

860  GOTO  890 

870  PLOT  K,LGT<T) 

880  NEXT  K 

890  PENUP 

900  FOR  K=Bins  TO  1  STEP  -1 

910  T=1.-T<K) 

920  IF  T>0.  THEN  940 

930  GOTO  960 

940  PLOT  K,LGT<T) 

950  NEXT  K 

960  PENUP 

970  PAUSE 

980  DUMP  GRAPHICS 

990  END 


SIGNAL 

NOISE 

NOISE 


SIGNAL 
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APPENDIX  E.  PROGRAM  FOR  COLORED  INPUT  NOISE 


One  program  is  listed  in  this  appendix.  It  computes  the 
covariance  matrices  and  eigenvalues  as  described  in  appendix  D, 
but  now  for  a  colored  noise  input.  The  program  for  computing  the 
receiver  operating  characteristics  from  the  eigenvalues  is 
identical  to  that  listed  above  and  therefore  has  not  been 
repeated  here. 

10  !  COMPUTE  COVflRIflMCE  MATRICES  AMD  EIGEHVALUES 

20  !  STORE  EIGEHVALUES  IH  "EIG"  IH  LIME  2260 


Ka=-10 
Kb  =  5 
Kc  =0 
Kd  =  5 

Del ta=. 2 
J  =  3 

DATA  11., -48. 
DATA  1 . , . 1 25 , 
M  =  2 

DATA  39. , 60. 
DATA  .2, .4 
N  =  4 

DATA  l.,2.,3 
DATA  . 3 , . 5 , . 


75. 

066667 


l.|2.,3.,4. 
.  3, . 5, .  7,  .  9 


INPUT  SIGNAL  START 
INPUT  SIGNAL  END;  Kb  >=  Ka 
ACCUMULATOR  START 
ACCUMULATOR  END;  Kd  >=  Kc,Ka 
TIME  SAMPLING  INCREMENT  c:  SECONDS; 
NUMBER  OF  SIGNAL  COMPONENTS 
SIGNAL  SCALINGS  <NE  SET  SUM  =  1) 
SIGNAL  TIME  CONSTANTS  (SECONDS) 
NUMBER  OF  NOISE  COMPONENTS 
NOISE  SCALINGS  (NE  SET  SUM  =  1) 
NOISE  TIME  CONSTANTS  (SECONDS) 
NUMBER  OF  F'LTER  COMPONENTS 
FILTER  SCALINGS  (ARBITRARY) 
FILTER  TIME  CONSTANTS  (SECONDS) 
NUMBER  OF  S I GNAL-TO-NO I SE  RATIOS 
,10  !  INPUT  SNRS  IN  DE 


Nr=20  !  NUMBER  OF  S I GNAL-TO-NO I SE  RATIOS 

DATA  0,1,2,3,4,5,6,7,8,9,10  !  INPUT  SNRS  IN  D1 

DATA  11,12,13,14,15,16,17,18,19 
IF  (Ka<=Kb)  AND  (Kc<=Kd>  AND  (Ka<=Kd>  THEN  230 
PRINT  "PROBLEM  WITH  PARAMETERS" 

STOP 

DIM  Al pha( 10) , A( 10) , Beta< 10) , B( 10) , Ps i ( 1 0 ) , C ( 1 0 ) , R ( 20 ) 
DIM  W( 100) , Pc ( 10) , Ec ( 10) ,Pcw< 10) , Eb( 10) , Chb( 10) , Shb( 10) 
DIM  Che  ( 10) ,  She  (  10) ,  Mu(  10)  ,  Gainina(  10) ,  Cn(200) ,  Ea(  10) 

DIM  E(0, 1000) , Gca( 100) , Gcda( 100) , Gcb( 100) , Gbda( 100) 

DIM  Cs(5000) , Cw(0, 10000),Ca(0, 1 0000 ), D( 1 00 ), E 1 g ( 0 , 2000 ) 
DOUBLE  Ka,Kb,Kc,Kd,J,M,N,Nr,Ns,Ks,Ms,Kdc,Ko 
DOUBLE  LI , L2, L, LI  1 , Js, Qs, Ki , Ks 1 , K 1 , KS2 , K2 , I 
REDIM  Alpha( 1 : J), A( 1 j  J) 

READ  Alpha(*), A(*) 

MAT  A 1 pha=A 1 pha^ ( SUM ( A 1 pha) ) 

MAT  A=(De1ta)/'R 
REDIM  Beta(liM),B(liM) 

READ  Beta(*),B(*) 

MAT  Bet a=Bet a/ ( SUM( Bet  a) ) 

MAT  B=(De1ta)/B 
REDIM  Psi ( 1 «  N) , C( 1  I N) 

READ  Psi (♦) , C(*) 

MAT  C=(Delta)/C 
REDIM  R( 1 : Nr) 

READ  R(*) 

REDIM  W(KciKd) 

CALL  Wei ghts<Kc , Kd, W(*) ) 

MAT  W*SQR(W) 
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460  REIUM  Pc  <  1  ;  M)  ,  Ec  <  1  :  M)  !  FILTER  OUTPUT  NOISE  COVFlF:  I  RHCE 

470  FOR  N£=1  TO  N 

480  C=C<N£) 

490  Pc  OIs)=P£i <N£)*C 

500  Ec <Ns)=EXP<-C) 

510  NEXT  Ns 

520  REDIM  PcMdiN) 

530  FOR  Ns=l  TO  N 

540  E=Ec<Ns) 

550  S=0. 

560  FOR  Ks=l  TO  N 

570  s=s+Pc(:ks) '<:i.-e*ec<:ks)> 

580  NEXT  Ks 

590  Pcv(Ns)=Pc  <:Ns)»S 

600  NEXT  Ns 

610  REDIM  Eb<l;M) 

620  FOR  M:=l  TO  M 

630  Eb(Ms)=EXP<-B<M5) ) 

640  NEXT  Ms 

650  REDIM  Chb< 1 : M) , Shb< 1 ! M) 

660  FOR  Ms=l  TO  M 

670  E=Eb<Ms) 

680  R=l./E 

690  Chb<Ms)=. 5*<R+£) 

700  Shb<:M£)  =  .  5*<R-E) 

710  NEXT  Ms 

720  REDIM  Che < 1 : N) , She <  1  j  N) 

730  FOR  Hs=l  TO  N 

740  E=Ec<Hs) 

750  R=l./E 

760  Che <Ns)=. 5*<R+E) 

770  She <Hs)=. 5*<R-E) 

780  NEXT  Ns 

790  REDIM  Mu<l:M) 

800  FOR  Ms=l  TO  M 

810  Chb=Chb<Ns) 

820  S=0. 

830  FOR  Ns=l  TO  N 

840  S  =  S  +  Pc'ANs)*Shc  C  Ns  )  /  <  Che  <  Ns  ) -Chb ) 

850  NEXT  Ns 

860  Mu < Ms > =Eet a< Ms ) *3 

870  NEXT  Ms 

880  REDIM  GamrnadzN) 

890  FOR  Ns=l  TO  N 

900  Che=Chc<Ns) 

910  S=0. 

920  FOR  Ms=l  TO  M 

930  S  =  S  +  Bet  a<Ns)*Shb<Ms)/'<Chb<Ms>-Chc  ) 

940  NEXT  Ms 

950  GammaC Ns ) =Pc u < Ns ) *S 

960  NEXT  Ns 

970  Kde=Kd-Kc 

980  REDIM  Cn-;0;Kdc) 

990  MAT  Cn=<0.) 

1000  FOR  Ms=l  TO  M 

1010  Nu=Nu(Ms) 

1020  E=Eb<M3) 

1030  P=l. 

1040  Crn(0)=Cn<0)+Mu 

1050  FOR  Ks=l  TO  Kdc 

1060  P=P*E 
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1070 
I  080 
I  090 
I  100 
1110 
1  1  20 
1  1  3  0 
1  140 
1  1  5  0 
1  160 
1  1  70 
I  ISO 
1  1  90 
1  200 
1210 
1220 
1  2  3  0 
1240 
1  250 
1260 
1  2  7  0 
1280 
12  9  0 
1  3  0  0 
1  3  1  0 
1320 
1  3  3  0 
1  34  0 
1  350 
!  360 
1  3  70 
1  380 
I  3  9  0 
1  400 
14  10 
14  20 
14  30 
1440 
1  4  o  0 
14  6  0 
14  70 

14  6  0 
1  4  9  0 
1  500 

15  10 

1  j  2  0 

1  • 
1  5  4  0 
1550 

1  j  6  0 

15  70 
156  0 


Cna<5>=Cr-i<Ks)+r'1u*P 

NEXT  Ks 

NEXT  Ns 

FOR  Ns=l  TO  N 

G  0  Stfii ni  i  C  N  s  ^ 

E  =  Ec  <Ns> 

P=1  , 

C  ri  <  0  '^  =  C  fi  <  0  >  +  G  a 
FOR  Ks=l  TO  Kdc 
P  =  P*E 

C n  c  K  s  )  =  C  r'l  (  K  s  )  +  G  a  *  P 
NEXT  Ks 
NEXT  Ns 

FOR  Js=l  TO  J  !  FILTER  OUTPUT 

H  (  J  s  .>  =  E  X  P  <  fl  <  J  s  )  ) 

NEXT  Js 

l:  u  =  MflXO(a,  Kc  ) 

REDIN  Eac:  1  :  N  )  ,  E  <  1  :  N  ,  Ko  s  Kd  ) 

FOP  Ns=l  TO  N 
E=Ec (Ns) 

Ea(Ns)=l . /E 

E(Ns , Ko)=EXP<-C(Ns)*Ko) 

FOR  Ks=Ko+l  TO  Kd 
E(Ns,Ks)=E(Ns,Ks-l)*E 
NEXT  Ks 
NEXT  Ns 

Ll=NIN(Ko,Kb)+l 
L2=NIN(Kd, Kb)+1 

RE  D I N  Gc  a<  L 1 ! L2  > , Gc  da( L 1 : L2 ) , Gc  b  <  L 1 : L2 
L=Kd-Ko+ 1 

REIHH  Cs  <  1  :  L*<L+ D/2> 

L 1 1 =L 1  +  1 
FOR  Js=l  TO  J 
H  1  =  H  1  fj  Fi  i  (  J  S  ) 

H  =-  fl  (  J  S  > 

H  1  =  H  »  H  -  1  . 

FOR  Ns=l  TO  N 
C=Ea(Ns) 

Ca=C*fl 
L  da  =  C  /Fl 
0  s  1  =  1  .  ~  L  a 
H  c  -  FI  -  C 

Hs  1 ♦Pc ( Ns  ) 

6  c  ii  - 1  .  i  '  K  a  -’'  Cal 

Gc  S'  L  1  ,)=Ca  'Ll  /  Cal 

Gc  .)  ac  LI  )  =Cd  a-L  1  ■'  <  1  .  -Cda) 

F 0 P  l's  =  Lll  TO  L2 
G  c  a  I  F  s  )  —  L  a  ♦  i-j  c  a  ( 1  s  “  1  ) 

G  c  d  a  ■'  F  s  -'  =  C  d  a  ♦  G  c  d  a  (  K  s  -  1  ) 

NEXT  Ks 

FjK  03=1  TO  N 

b  ■-  E  a  <  Us) 


iIGNfiL  COVRRIhNCE 


,  G  b  d  a  (  L  1  :  L  2  ) 
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1680  Sl=Cb''Ka/Cbl*<Fll+Cbl  )^<Flb*flc> 

1690  GcbCLl  )=Cb"'Ll/Cbl 

1700  GbdaCLl )=Bda-Ll/< 1 . -Bda) 

1710  FOR  Ks=Lll  TO  L2 

1720  Gcb<Ks)=Cb*Gcb<Ks-l ) 

1730  Gbda< Ks ) =Bda*Gbda< Ks- 1  > 

1740  NEXT  Ks 

1750  Ki=0 

1760  FOR  Ksl=Ko  TO  Kd 

1770  Kl=MIN<Ksl , Kb)+1 

1780  S2=Sl+Gcb(Kl)*F-Gba*Gcda<Kl > 


1790 

S3  —  GcaCKl )— Gca 

1800 

E  =  Aee*ECN£, Ksl  ) 

1810 

FOR  Ks2=K£l  TO  Kd 

1820 

Ki =K1 +1 

1830 

K2=MINCKs2, Kb)+1 

1840 

S=S2+GbdaCK2)*S3  ! 

SUM  S 

1850 

CsCK'i  )=C£CKi  )+E*ECQ£,Ks2)*S 

I860 

NEXT  KS2 

1870 

NEXT  Ksl 

1  880 

NEXT  Qs 

1890 

NEXT  Ns 

1900 

NEXT  Js 

1910 

L  =  Kd-K:c  +  l  ! 

TOTAL 

WEIGHTED  COVARIANCE  MATRIX 

1920 

REDIM  CwC 1 : L, 1 : L) , CaC 1 

:L, 1:L), 

DC  1 : L ) , E 1 g  C  0 : Nr , 1 ; L ) 

1930 

LI l=Kc-l 

1  940 

R  =  0. 

1  950 

FOR  1=0  TO  Nr 

1  960 

K  1  =  0 

1970 

IF  1=0  THEN  1990 

1  980 

R=10. -'CRC  I  )*.  1  )  ! 

INPUT 

POWER  SIGNAL-TO-NOISE  PATIO 

1990 

FOR  K£l=Kc  TO  Kd 

2000 

W=WCKsl)  ! 

WEIGHTS  <wCk)> 

2010 

FOR  Ks2=Ksl  TO  Kd 

2020 

IF  K5l>=Ko  THEN  2050 

2030  C£=0. 

2040  GOTO  2070 

2050  K1=K1+1 

2060  Cs=CsCKi) 

2070  Pr  =  W*N<Ks2^*'.';Cn<K32-Ksl  >+R*Cs) 

2080  Ll=Ksl-Lll 

2090  L2=Ks2-Lll 

2100  Cw(Ll ,L2)=Cw(L2,Ll)=Pr 

2110  NEXT  Ks2 

2  120  NEXT  k'sl 

2130  CALL  SudC L , L , ♦ ) , C *) , DC *> )  '  EIGENVALUES 

2140  MAT  SORT  DC*)  DES 

2150  IF  DCL>>0,  THEN  2180 

2160  PRINT  "PROBLEM:  SOME  NON-POSITIVE  EIGENVALUES" 

2170  PAUSE 

2180  PRINT  "I  ="!l;"  CONDITION  NUMBER  ="(DclV'D<L> 

2190  PRINT  DC*) 

2200  PRINT 

2210  FOR  K£=1  to  L 

2220  E  1  gC  I  ,  K£)=D'^K£  )  !  STORE  EIGENVALUES 

2230  NEXT  Ks 

2240  NEXT  T 

2250  NASS  STORAGE  IS  ":CS80,7" 

2260  ASSIGN  #1  TO  "EIG" 

2270  PRINT  # 1 ; Nr , Kc , KU, E 1 gC ♦  ) 

2  280  A  S  S  i  G  r  i  1  T  u  * 

2290  END 
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